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ABSTRACT 


Circular-to-rectangular transition ducts are used on aircraft with rectangular ex- 
haust nozzles to connect the engine and nozzle. Often, the incoming flow of these 
transition ducts includes a swirling velocity component remaining from the gas tur- 
bine engine. Previous transition duct studies have either not included inlet swirl or 
when inlet swirl was considered, only overall performance parameters were evaluated. 
This dissertation reports the study of circular-to-rectangular transition duct flows with 
and without inlet swirl. This study was completed in order to understand the effect 
of inlet swirl on the transition duct flow field and to provide detailed duct flow data 
for comparison with numerical code predictions. 

A method was devised to create a swirling, solid body rotational flow with min- 
imal associated disturbances. Details of the swirl generator design and construction 
are discussed. Coefficients based on velocities and total and static pressures mea- 
sured in cross stream planes at four axial locations within the transition duct along 
with surface static pressures and surface oil film visualization are presented for both 
nonswirling and swirling incoming flows. For nonswirling flow, measurements were 
recorded for three inlet centerline Mach numbers; 0.20, 0.35, and 0.50. The cen- 
terline Mach number for the swirling flow measurements was 0.35. The maximum 
ideal swirl angle was 15.6°. 

A method was developed to acquire trace gas measurements within the transition 
duct at high flow velocities. Trace gas results are presented for the case of nonswirling 
incoming flow with an inlet centerline Mach number of 0.50. Statistical methods are 
used to help interpret the trace gas results. 

Inlet swirl significantly changes the transition duct flow field. For nonswirling 
flow, the distribution of static pressure is attributed to the response of the flow field 
outside the boundary layer to the changing duct geometry. With swirl, the static 
pressure distribution is additionally affected by the streamline curvature of the swirl. 
For nonswirling incoming flow the static pressure distribution produced skew-induced 
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secondary flows within the boundary layer that evolved into two pair of counter- 
rotating side wall vortices. These vortices were absent in the swirling incoming flow 
case. The results show the effects of inlet swirl should be included in the design 
and numerical analysis of circular-to-rectangular transition ducts for aircraft exhaust 
systems. 
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CHAPTER I. INTRODUCTION 

Ducts are used in aircraft propulsion systems to conjoin different components. 
For example, S-shaped ducts are used to connect the inlet and engine of Boeing 727 
and LotMieedi^lOM’airei^.^ ^Recently, nonaxisymmetric exhaust nozzles have been 
proposed for aircraft applications which offer improved maneuverability by means of 
thrust vectoring and reduced thermal emission in comparison to axisymmetric nozzles. 
Rectangular exhaust nozzles are visible on the F-117A fighter and the B-2 bomber* 
prototype. Aircraft with rectangular exhaust nozzles require a circular-to-rectangular 
transition duct to connect the engine and nozzle. 

Detailed aerodynamic measurements of the flow field of a circular-to-rectangular 
transition duct are important for at least two reasons; to identify the important fluid 
dynamic phenomena which determine the transition duct flow field, and to provide 
data to compare with CFD predictions of the transition duct flow field. Knowledge of 
the relevant fluid dynamics can guide the engineer in improving the design of circular- 
to-rectangular transition ducts. Comparison with experimental data can provide a 
benchmark to measure the accuracy of CFD flow field predictions which lead to 
improvements in CFD methods. 

To the designer, the value of experimental model data is enhanced when experi- 
mental test conditions accurately reproduce the conditions encountered in applications. 
For circular-to-rectangular transition ducts the incoming flow often includes a signif- 
icant rotational velocity component, which remains from the engine turbine. Inlet 
s wirl can significantly .alter rhe flaw field throughout the transition duct. 

This dissertation reports the study of circular-to-rectangular transition duct flows 
with and without inlet swirl. The objective of this research was to ascertain the effect 
of inlet swirl on the transition duct flow field and to provide detailed duct flow data 
for comparison with CFD flow field predictions. Previous transition duct studies have 
either not included inlet swirl, or when swirl was considered only overall performance 
parameters were evaluated. 
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An improved understanding of circular-to-rectangular transition duct fluid dy- 
namics was gained though detailed measurements of the transition duct flow field. 
Coefficients based on extensive measurements of velocity, total pressure and static 
pressure, acquired near the inlet and in four cross stream measurement planes within 
a circular-to-rectangular transition duct, are presented for flows with and without inlet 
swirl. Surface static pressure coefficients and surface oil film visualization results are 
also displayed. Additionally, a trace gas technique was used to obtain a Lagrangian 
description of the transition flow field without inlet swirl and those results are reported. 

The simulation of circular-to-rectangular transition duct flows was improved 
through the use of a new swirl generator conceived, designed and constructed for 
this study. This swirl generator provides a solid body rotational flow with minimal 
disturbances, suitable for detailed experimental measurements and comparison with 
CFD predictions. A novel and original method was invented to produced solid 
body rotational flows which employs blades fabricated from cylindrical tubing with 
diminishing chord length, followed by rotating honeycomb and screen. The design, 
testing and operation of this swirl generator are reported herein. 

A refinement of the five-hole probe aerodynamic measurement calibration and 
data reduction technique was used extensively for transition duct flow field measure- 
ments. This measurement technique is based on an analysis that provides a mathe- 
matical foundation for the calibration and use of five-hole probes. The results of the 
mathematical analysis used to develop the calibration and data reduction procedures 
and the numerical algorithms used to implement these procedures are included in an 
appendix. 

The ethylene trace gas technique was adapted to higher speed (e.g., compressible, 
subsonic) flows and used to obtain Lagrangian data in the circular-to-rectangular 
transition duct. Trace gas results are presented for the flow field without inlet swirl. 
Statistical methods are used to quantitatively interpret the trace gas results. Numerical 
algorithms that implement the procedures used to calculate trace gas statistics are 
included in an appendix. Novel methods are described that were employed to inject 



ethylene trace gas in the transition duct flow field with minimal disturbances and 
to sample the flow with effectively unlimited spatial resolution of the measurement 
locations within a cross stream plane. 

The experimental measurements and analysis of the circular-to-rectangular tran- 
sition duct flow field reported in tins dissertation will be valuable to an engineer 
designing circular-to-rectangular transition ducts or using CFD methods to calculate 
similar internal flows. The methods described and used for swirl generation, five- 
hole probe calibration and data reduction, and ethylene trace gas measurements will 
be useful to an engineer employed in similar experimental studies. 
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CHAPTER II. REVIEW OF PRIOR RESEARCH 

Past work is briefly reviewed in this chapter to offer a perspective on how the 
present research project was conceived and why it was conducted as reported. This 
project is an extension of progress made in past development of transition ducts. 
Transition duct research is discussed first and the case for new data is established. 
The need for an effective swirl generator is ascertained next. Finally, the usefulness 
of the trace gas technique as a vital measurement tool in this project is proposed and 
the adaptation of this technique to compressible flow is described. 

Transition Duct Research 

Nonaxisymmetric exhaust nozzles are employed on high performance aircraft to 
improve performance. The Pratt and Whitney two-dimensional/convergent-divergent 
nozzle described by Stevens et al. [1] is a rectangular nozzle designed to allow 
thrust vectoring, thrust reversal, jet area variation, and expansion area variation. 
Reduced thermal emission is another benefit attributed to rectangular exhaust nozzles. 
Rectangular exhaust nozzles require a circular-to-rectangular transition duct to connect 
the engine and nozzle. To allow proper exhaust nozzle performance, the transition duct 
should minimize transverse velocity components at the duct exit and total pressure 
losses throughout the duct. 

The experimental study of circular-to-rectangular transition duct flows is one 
component of several in an effort at the NASA Lewis Research Center to accumulate 
experimental measurements of subsonic flows through various duct geometries. These 
data are used to identify and understand primary fluid phenomena governing the flow 
fields involved and to compare with relevant CFD predictions. These experimental 
efforts support the extensive enterprise of developing CFD methods to aid the design 
and analysis of aircraft propulsion components and systems. A review of the related 
numerical methods development is given by Anderson [2]. 
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Under grants awarded by the NASA Lewis Research Center, Taylor et al. acquired 
laser Doppler velocimetiy measurements in a square cross section duct with 90° bend 
[3], a square-to-round transition duct [4], an S-shaped square cross section duct [5] 
and an S-shaped circular cross section duct [6]. Enayet et al. [7] made comparable 
measurements in a circular cross section duct with a 90° bend. A common feature 
of these five studies is that measurements were made for both laminar and turbulent 
incompressible flow with thin inlet boundary layers. Also funded as part of the 
NASA Lewis Research Center effort were studies by Vakili et al. of compressible, 
flow through an S-shaped circular cross section duct with constant cross-sectional 
area [8-10] and compressible flow through an S-shaped circular cross section duct 
with increasing cross-sectional area[Tl-T3]. The research of Vakili et al. was also 
characterized by thin inlet boundary layers. 

Researchers have also experimentally explored the aerodynamics of circular-to- 
rectangular transition ducts. Patrick and McCormick [14, 15] acquired experimental 
measurements at the inlet and exit planes of two circular-to-rectangular transition 
ducts, the first with an aspect ratio of 3 and a length-to-diameter ratio of 1 and the 
second with an aspect ratio of 6 and a length-to-diameter ratio of 3. The first duct 
maintained a constant cross-sectional area throughout the transition Tegion. In the 
second duct the cross-sectional area at the inlet and exit of the transition region were 
equal but the cross-sectional area increased up to 10% in the transition region. Total 
pressures were measured with pressure probes and the three mean velocity and three 
normal Reynolds stress components were measured with a laser Doppler velocimeter. 
Eor this research the inlet Mach number was nominally 0.10, the Reynolds number, 
based on the inlet velocity and the duct inlet diameter, was 460,000 and 320,000 for 
the first and second ducts, respectively. The difference in Reynolds number resulted 
from the ducts having different inlet diameters. The ratio of the inlet boundary layer 
thickness to the duct inlet diameter was approximately 0.10 for both ducts. 

Miau et al. [16, 17] acquired measurements in three circular-to-rectangular 
transition ducts for three different inlet flow conditions. The aspect ratio of each 
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duct was 2.0, the length-to-diameter ratios of the transition region were 1.08, 0.92, 
and 0.54. The inlet Mach numbers ranged between approximately 0.018 for the 
lowest inlet velocity condition to 0.106 for the highest inlet velocity. The Reynolds 
numbers, based on the free stream velocity and the duct inlet diameter, ranged between 
195,000 and 1,140,000. The ratios of the inlet boundary layer thickness to the duct 
inlet diameter ranged between 0.0205 and 0.0175. Mean velocities at the transition 
duct exits and within the boundary layer along the diverging walls and turbulence 
intensity at the transition duct exits were measured with a hot wire anemometer. 
Surface static pressures were measured along 12 cross stream planes and water tunnel 
dye visualization results were also presented. 

In a benchmark study, Davis and Gessner [18] measured static pressures, mean 
velocities, and all six Reynolds stress components in three cross stream planes within 
a circular-to-rectangular transition duct. Mean velocity and turbulence measurements 
were made with a hot wire anemometer, static pressures were measured with a 
pressure probe. The ratio of the boundary layer thickness to the inlet duct diameter 
was approximately 0.14. The Reynolds number, based on the bulk velocity and duct 
inlet diameter, was 390,000. The free stream Mach number was nominally 0.10. The 
transition duct of Davis and Gessner was identical to the one used in this study, and 
is described in detail in Chapter in. 

In exhaust component applications the incoming flow to the circular-to-rectangular 
transition duct is turbulent and subsonic. Often, the incoming flow retains a significant 
rotational velocity component. The term swirl refers to this rotational velocity 
component, which remains from the engine turbine. Representative studies of turbine 
exit flow angles [19-21] have shown that swirl often exists at turbine design operating 
conditions, and may be as great as 30° or more at off-design operating conditions. 
Inlet swirl can significantly alter the flow field throughout the transition duct. 

Each of the above mentioned circular-to-rectangular transition duct studies in- 
volved incompressible flow without inlet swirl. Burley et al. [22, 23] tested five 
circular-to-rectangular transition ducts including one with swirl vanes installed up- 
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stream of the duct inlet. Four of these transition ducts maintained a constant cross- 
sectional area throughout the transition region. For one transition duct the cross- 
sectional area decreased in the downstream direction. The aspect ratio of all five 
transition ducts was 6.33. The measurements were, however, limited to surface static 
pressures, a thrust ratio performance parameter and a discharge coefficient. No infor- 
mation was provided regarding the inlet Mach number, inlet boundary layer thickness 
or Reynolds number for these experiments. Test conditions were described by a noz- 
zle pressure ratio. 

The objective of the transition duct research reported in this dissertation was 
dtwo-fold; to ascertain the ^effect of inletmdrliontire tiai^ 
provide a set of detailed data for comparison with CFD predictions of transition duct 
flows with and without inlet swirl. A swirl generator was designed and constructed 
to add a swirling component to the flow immediately upstream of the transition duct. 
The intent of the swirl generator was to create a solid body rotational flow, free of 
wakes and other disturbances. 

Swirling Flow Generation 

Ahmed et al. [24] studied the generation and management of swirling flows. The 
focus of their effort was to study how flow manipulators (e.g., honeycomb, screen, 
etc.) affect swirling flows. Three different swirl generator designs were investigated. 
The first, called the airfoil swirl generator, consisted of two NACA 0012 airfoils 
mounted in a circular cross section duct. Each airfoil was approximately one pipe 
Tadius long, allowing the two airfoils to join at a common hinge point at the centerline 
and span the duct. Swirl was induced by inclining the airfoils to the centerline. This 
swirl generator produced a strong rotational core. Plots of the radial distribution of 
axial vorticity, measured with a vane vorticity indicator, showed the maximum axial 
vorticity occurred at the flow center and reached zero near the wall. A tangential-jet 
swirl generator was the second configuration tested. This apparatus produced swirl 
by injecting air through eight equally spaced tangential jets into the primary stream 



8 


of a circular cross section duct. This configuration produced strong rotation near 
the walls. Measurements showed the radial distribution of axial vorticity reached a 
maximum very near the wall and rapidly approached zero a short distance away from 
the wall. The final swirl generator was the named the swirling-jet ejector. Swirl 
was generated by flow through a section of annular vanes followed by a radial flow 
section. The swirling flow was used as the primary flow of an ejector. Like the 
airfoil swirl generator, this device produced a strong rotational core within a nearly 
irrotational outer region. 

The swirl generator used by Burley et al. [22, 23] in their transition duct studies 
employed 12 flat plate vanes attached to a central hub and inclined to the flow axis. 
This swirl generator is similar in principal to the airfoil swirl generator of Ahmed et 
al [24], although no measurements of the flow field exiting the swirl generator were 
presented. The swirling flow cases studied were characterized by the angle the vanes 
were inclined to the centerline. 

Hallett [25] described a swirl generator designed for the investigation of jet 
mixing and vortex breakdown. Air was injected tangentially from four slots along 
the periphery of the principal stream. This swirl generator is similar in principal to 
the tangential-jet swirl generator of Ahmed et al. 

A rotating pipe was used to generate swirling flows by Weske and Sturov [26]. 
The purpose of their research was to study the turbulence of confined swirling 
flows. Their measurements show this swirl generator produced a flow with solid 
body rotation. The power required to rotate the pipe was provided externally. The 
facility size and mass flow rates of air through the facility were low compared to the 
requirements for aircraft transition duct studies. 

The review of swirl generation literature revealed a new swirl generator design 
was necessary for the transition duct study. The swirl generator designs studied by 
Ahmed et al., and their variations used by Burley et al. and Hallett do not produce the 
desired solid body rotation in the flow. The swirl generator of Weske and Sturov is 
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promising, it produces a flow with solid body rotation but requires an external motor 
to supply the power to rotate the pipe. A new design which employs the passive 
features of the airfoil type swirl generator with a rotating pipe was developed. The 
design of the new swirl generator is described in Chapter HI, and measurements of 
the swirling flow field are presented in Chapter V. 

The Trace Gas Technique 

The trace gas technique involves injecting and tracking a nonreacting discernible 
foreign gas in the host flow field being investigated. Information about the host 
flow field is acquired by measuring the concentration of the trace gas in the host 
gas at many locations downstream of the point of injection. The trace gas technique 
is a valuable tool for several reasons; it is relatively inexpensive, easy to use, and 
is a source of Lagrangian flow field data. Previous trace gas experiments can be 
characterized by several important features, such as: 

1. The objectives the trace gas technique were used to achieve 

2. The types of trace gas used and their methods of detection 

3. The class of flow fields investigated 

4. The methods of trace gas injection used 

5. The trace gas sampling methods employed 

6. Important information about the trace gas technique that were revealed by the 
experiments 

The relevant trace gas literature will be. discussed in the context of these features. 

Previous experiments can be divided into two broad categories distinguished by 
the objectives the trace gas technique were used to achieve. The two classifications 
are inferential experiments and direct experiments. Inferential experiments refer 
to studies where the information about the flow field was deduced from the trace 
gas concentration measurements. Inferential trace gas experiments often require a 
priori assumptions about the host flow field being investigated or additional flow 
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field data acquired by other measurement techniques. The information deduced from 
the concentration measurements may be the mean streamline of a fluid particle or 
the relative importance of convection and diffusion in the mixing of trace gas with 
the host flow. 

Some important inferential trace gas experiments include the research of Kerre- 
brock and Mikolajczak [27], who tracked fluid flow from the pressure side of a rotor 
blade through the downstream stator row. Denton and Usui [28] studied the transport 
of fluid from the end wall boundary layers of a single-stage low speed air turbine. 
Smith [29], Moore [30], and Moore and Smith [31] used the trace gas technique in 
a linear turbine cascade to help locate the source of high and low total pressure loss 
fluid. Gallimore [32] and Gallimore and Cumpsty [33] used the trace gas technique 
in two low speed, four-stage air compressors. Their investigation was performed to 
ascertain the relative importance of convection and turbulent diffusion in spanwise 
(radial) mixing within multistage compressors. Bauer [34] and Wisler et al. [35] also 
used the trace gas technique to investigate spanwise convection and mixing in the 
third stage stator row of a low speed, four-stage, axial flow research compressor. 

Direct experiments refer to studies where the extent of mixing of the trace 
gas into the flow field was of interest. In this case the final result was made 
by direct observation. Often, the trace gas was used to represent another gas, or 
the trace gas concentration represented another fluid property such as temperature. 
Direct experiments include the work reported by Hallett and Guenther [36], who 
used the trace gas technique to quantify the degree of mixing within a combustor. 
McGreehan et al. [37] employed the trace gas technique to measure the rate at which 
airflow leaked through a labyrinth seal in a gas turbine engine. Joslyn and Dring 
[38] performed an experiment in a large scale low speed turbine where trace gas 
concentration was used as an analogy to temperature to study the radial redistribution 
of temperature in the turbine. Folayan and Whitelaw [39] and Humphrey and 
Whitelaw [40] used the trace gas and temperature analogy to better understand film 
cooling methods on convex and concave surfaces. 
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Ethylene was used for the trace gas in References 28-35. Other gases that 
have been used as a trace gas are helium [27, 37, 39, 40], carbon dioxide [37, 38], 
and natural gas [36]. The instruments used for trace gas detection were largely 
determined by the choice of trace gas. A flame ionization detector provides a very 
accurate measurement of ethylene concentration and was used in References 28-35. 
Researchers that use helium have frequently used a thermal conductivity cell for trace 
gas detection [27, 37, 39, 40]. Carbon dioxide concentration was commonly measured 
with a nondispersive infrared detector [37, 38]. 

The class of flows investigated by the trace gas technique may be categorized 
by the geometry or application Gf the experiment &uch as iurbomachines or mixing 
nozzles, or more importandy by features of the flow, particularly the flow velocity. 
All of the uses of ethylene trace gas that have been cited above have been at relatively 
low velocities, in nearly all cases the Mach number of the host flow was less than 
0.10. One objective of the research reported in this dissertations was to adapt the 
ethylene trace gas technique to higher speed, compressible, subsonic flows. Building 
on this success, Davis et al. [41] recently used the ethylene trace gas technique at 
NASA Lewis Research Center to measure the mixing of two supersonic streams in 
a study of mixing nozzles for scramjet applications. 

The methods used to inject the trace gas into the host flow field can be classified 
as either surface injection methods or probe injection methods. The choice of trace 
gas injection method largely depends on the objectives of the trace gas experiment. 
Static pressure taps were used for surface trace gas injection in References 27, 28, 34, 
35. Slot shaped openings have also been used for surface trace gas injection [39, 40]. 
The most common trace gas probe injection method used was an L-shaped probe, like 
a downstream facing total pressure probe. This type of probe, or variations of it, have 
been used in References 28-35. The diameters of the injection probes reported in 
these references range from 3.2 mm to 0.56 mm. The results of experiments designed 
to measure and minimize disturbances from L-shaped type probes are reported in 
References 29, 30, 32, 33 and by Moore et al. [42]. These results are discussed 
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in Chapter IV. A different approach to trace gas injection was employed by Joslyn 
and Dring [38]. In their experiment trace gas was injected simultaneously at many 
point locations distributed throughout the injection plane to simulate the temperature 
distribution in an air turbine. 

The most common trace gas sampling probe design used in the past employed 
features similar to a single-hole Pitot type pressure probe [27-31, 34, 35]. The flow 
field sample was drawn through the single probe opening. Another trace gas sampling 
probe design reported was similar to a three-hole “cobra” type pressure probe [32, 
33]. When using this kind of probe, the two side holes were used to align the yaw axis 
of the probe with the flow field direction while the center hole was used to draw the 
sample. Other sampling methods involved the use of a trace gas sampling probe rake 
[29-31, 38] and sampling through surface static pressure tape openings [39, 40, 38]. 

Several ethylene trace gas experiments have revealed important information about 
the ethylene trace gas technique. The paper of Wisler et al. [35] generated consid- 
erable discussion. Some controversy centered around their conclusions regarding the 
relative importance of convection and turbulent diffusion in the radial transport of 
fluid in multistage axial flow compressors. Differences in the interpretation of the 
trace gas results, that is, deciding the roles of convection and turbulent diffusion in 
producing the measured ethylene concentration distributions were proposed by discus- 
sants. For example, Gallimore and Cumpsty proposed different interpretations of the 
Wisler et al. data. Leylek and Wisler [43] continued the study reported by Wisler et 
al. Additional ethylene trace gas data and Navier-Stokes code computational results 
were used to aid trace gas data interpretation. The difficulties of interpreting trace 
gas concentration measurements are discussed in Chapter VI. 
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CHAPTER III. FACILITIES 

This chapter describes the experimental facilities used during the course of this 
research. The calibration tunnel and the swirl generator are described in detail in the 
first and fourth sections of this chapter, titled Calibration Tunnel and Swirl Generator, 
respectively. These facilities were designed and constructed specifically for this 
research. The Internal Fluid Mechanics Facility, described in the second section 
of this chapter, and the transition duct, described in the third section, already existed 
prior to the beginning this research. 


Calibration Thnnel 

A small air flow tunnel was designed and constructed at the NASA Lewis 
Research Center to test the trace gas equipment and measurement technique before 
trace gas measurements were attempted in the transition duct. This tunnel was 
designed to operate at flow velocities comparable to those used in the transition 
duct. The tunnel involved an air jet stream that expands into a quiescent chamber. 
The jet was designed to have a large core of uniform, mean, axially directed velocity 
with low turbulence intensity. The use of this tunnel and the results of the trace gas 
experiments conducted in this tunnel are described in Chapter V. 

Components 

The tunnel, pictured in Figure III. 1 consists of three principal components. They 
are: 

1. The setding chamber 

2. The contraction 

3. The test section 




Figure III. 1 The calibration tunnel used to test the 
trace gas equipment and measurement technique 

Settling chamber The settling chamber consisted of a 91.4 cm length of 25.4 
cm diameter schedule 40 pipe. The settling chamber employed a honeycomb and 
screen combination to reduce the free stream turbulence. Scheiman [44] observed 
that a honeycomb alone is more effective at reducing lateral turbulence while a screen 
alone is more effective at reducing axial turbulence. Scheiman concluded that the 
combination of screens downstream of honeycomb is a productive way to utilize the 
turbulence reduction characteristics of each device. 

Loehrke and Nagib [45] attribute the effectiveness of a honeycomb at reducing 
lateral turbulence to the constraint of the lateral velocity components by the honey- 
comb side walls. Accordingly, they recommended that the Reynolds number, based 
on the honeycomb cell diameter, be less than the value for transition to turbulent 
flow. The honeycomb used in the settling chamber was constructed of stainless steel 
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and had a 0.95 cm cell diameter. When the jet is operating at a choked condition 
(its maximum velocity), the air velocity within the settling chamber is approximately 
4.45 m/s, resulting in a Reynolds number of approximately 2,500. At all operating 
conditions but the highest velocities this Reynolds number is less than the value of 
Reynolds number generally accepted as indicating transition to turbulent flow, 2300. 

The experiments of Loehrke and Nagib also indicate the ratio of the honeycomb 
cell length to the cell diameter has a significant effect on the ability of the honey- 
comb to reduce turbulence. Of the three different lengths they tested, the shortest, 
with a cell length to diameter ratio of eight, was most effective at reducing turbu- 
lence. The honeycomb used in the calibration tunnel has a cell length of 5.08 cm, 
providing a length-to-diameter ratio of 5.3 which conforms with Loehrke and Nagib’s 
recommendation of using a relatively short length honeycomb. 

Schubauer et al. [46] showed that a wire screen is most effective in reducing 
turbulence when the Reynolds number, based on the screen wire diameter, is below 
a critical value. Reynolds numbers above this value produce vortex shedding down- 
stream of the screen, resulting in higher, rather than lower, turbulence intensities. This 
critical Reynolds number depends upon the screen solidity a , which is the ratio of the 
closed area of the screen to the total area. An approximate relation, based on data in 
Reference 46 is Ke cr j t = 77.5 -62.5a for solidity values in the range 0.2 < a < 0.5. 

The screen used in the calibration tunnel was a 24 mesh stainless steel wire 
screen with a solidity of a — 0.327 and a wire diameter of 0.1905 mm. Based on 
this solidity and wire diameter, the critical Reynolds number should be approximately 
57. When the jet is operating at a choked condition, the resulting Reynolds number 
is approximately 52. Therefore, in all operating conditions the Reynolds number is 
less than the critical value. 

Contraction An axisymmetric contraction with a 52:1 area ratio was located 
immediately downstream of the settling chamber. Contractions are used downstream 
of settling chambers to increase the mean velocity, lower turbulence intensities, and 
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insure uniform mean flow. This contraction was machined from a solid aluminum 
round bar. The shape of the inside radius of the contraction was machined to match 
the curve of the two smoothly joined cubics given by Equation in.L This curve is 
smooth in the sense it is continuous and has a continuous first derivative. In a study 
of axisymmetric wind tunnel contractions. Morel 147] compared six different families 
of curves and concluded that two smoothly joined cubics was the curve best suited to 
use as the wall shape for an axisymmetric wind tunnel contraction. In Equation III. 1 
R\ is the upstream radius of the contraction, R 2 is the downstream radius, and X is 
the length of the contraction. The values used for the construction of the contraction 
were &i = 12*541 son, ifo — 1.7145 cm and X = 25.083 cm. 


_ L= { 1 - 4 (i-f)(f) 3 °<i<5 

Rl 1 ft “ 4 ( 1 " ft)(i “ 1 )' 1 5 < f < 1 


(ni-i) 


The outside diameter of the contraction was a constant 25.4 cm, to allow the 
contraction to be installed inside a 30.5 cm length of 25.4 cm diameter schedule 40 
pipe. This section of pipe was drilled and tapped to accept set screws to secure the 
contraction section. The outlet of the contraction was connected to a short section of 
3.81 cm diameter stainless steel tube which exhausted as a jet into the test section. 


Test section The test section was constructed from a 25.4 cm diameter schedule 
40 pipe cross on which two additional legs of 25.4 cm diameter schedule 40 pipe were 
added, providing six legs. The contraction section was attached to one leg, while the 
opposite leg was connected to an exhaust line. Windows, constructed of aluminum 
frames and 15.2 cm diameter quartz lenses, were installed in opposite legs to provide 
optical access of the test section. The optical access allowed visual inspection of the 
trace gas injection and sampling probes under operational conditions. The remaining 
pair of opposite legs were blanked off with aluminum plates and were used to allow 
instrument access. 
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Internal Fluid Mechanics Facility 

All transition duct experiments were conducted at the NASA Lewis Research 
Center using the Internal Fluid Mechanics Facility. This facility, shown in Figure 
III.2, was designed to support the research of a variety of internal flow configurations. 
Air may be supplied from either the test cell or from a continuous source of pressurized 
air. After passing through a large settling chamber containing a honeycomb and two 
screens the flow proceeds through an axisymmetric contraction having an area ratio 
of 57:1. The flow passes from the contraction through either a straight pipe to provide 
a nonswirling uniform incoming flow to the transition duct or the swirl generator to 
provide a swirling incoming flow to the transition duct. After passing through the 
transition duct the flow is exhausted into a discharge plenum which is continuously 
evacuated by central compressor facilities. In Figure III.2 the Internal Fluid Mechanics 
Facility is shown with the transition duct installed, the straight pipe in place instead 
of the swirl generator, and with air drawn from the test cell through the open end 
of the settling chamber. 

Throughout this dissertation, discussion of the swirl generator or the straight pipe, 
which are upstream of the transition duct, are with respect to a cylindrical coordinate 
system, x,$,r. The rr-axis of the coordinate system is coincident with the swirl 
generator or straight pipe centerline. Distances have been nondimensionalized with 
respect to the swirl generator or straight pipe radius, 10.21 cm. The radius, rather than 
the diameter, was used to conform with standard convention for axisymmetric flows. 

For all transition duct experiments the air was supplied from the test cell. This 
configuration does not allow the Mach number and Reynolds number to be varied 
independently, their relationship is given in Equation III.2. 


Re = mR M (i+i^MA 


-r-7 


( 111 * 2 ) 


HO \ 2 

The total pressure and temperature of the airflow is dependent on the ambient pressure 
and temperature within the test cell. Figure HI. 3 is a plot of Equation III.2 for 
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Figure in. 3 The Mach number - Reynolds number 
relationship for the Internal Fluid Mechanics Facility 

conditions evaluated at the inlet of the transition duct. The example test cell air 
pressure and temperature for generating Figure III.3 are 101.3 KPa and 294 °K. The 
Reynolds number, shown on the ordinate of Figure III. 3 is per meter of characteristic 
length. The solid line represents nonswirling flow, the broken line represents swirling 
flow. The difference between these two curves results from total pressure losses 
produced by the swirl generator. These total pressure losses are discussed in the 
Swirl Generator Validation section of Chapter V. 

Transition Duct 

Figure in.4 shows a wireframe drawing of the circular-to-rectangular transition 
duct that was tested. All discussion of the transition duct throughout the dissertation 
is with respect to the Cartesian coordinate system shown in Figure III. 4. The origin 
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Figure III.4 Circular-to-rectangular transition duct 

of the coordinate system is located one diameter upstream of the beginning of the 
transition region. This origin coincides with the physical beginning of the transition 
duct (or the exit of the straight pipe or swirl generator). The z-axis is coincident to 
the duct centerline. Distances have been nondimensionalized with respect to the inlet 
diameter of the transition duct, 20.42 cm. 

This transition duct is a member of a family of ducts designed at the NASA 
Lewis Research Center. In the yz-plane through each cross section, the surface of 
the family of ducts satisfies Equation III.3. 



The parameters yo, r y , n y , zq, r z , n z , which specify the exact geometry of the 
members of this family, are all functions of the axial distance x . 


F : x -¥ yo,r y ,n y ,zo,r z ,n z 


(ni.4) 
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The functional relation represented by Equation III.4 is summarized in Table IILl. 

The transition duct geometry possesses two planes of symmetry, the xy-plane 
and xz-plane. If the point ( x,y,z ) satisfies Equation HI.4, then so does the points 

(x,-y,z), (x,y, —z), and ( x,-y,-z ). 

The transition region of the duct is the region where the cross section in the yz- 
plane varies with axial distance. The transition region lies within 1.0 < x/D < 2.5. 
The ratio of the length of the transition region to the inlet diameter is L/D = 1.5. 
The ratio of major to minor axis dimensions at the duct exit, referred to as the- 
aspect ratio, is r y /r z = 3.0. The cross-sectional areas at the inlet and exit are 
equal. In the transition region the cross-sectional area increases to 1.15 times the 
inlet area. The distribution of cross-sectional area in the transition region is equal 
to the cross-sectional area that would result if a duct with the same aspect ratio and 
transition region length to inlet diameter ratio were constructed of circular and flat 
sections of metal only. This type of construction was believed to more accurately 
represent fabrication methods used in the aviation industry. The cross-sectional shape 
of these transition ducts becomes more rectangular as the 'value of the exponents 
increase, but they are never truly rectangular. This was done to provide the smooth 
boundary required by some CFD methods. The circular-to-rectangular ducts studied 
in References 14, 15, 18, 22, 23 are all members of this family of transition ducts. 

The transition duct was constructed in two halves joined at the xy-plane. Each 
half was machined from a single piece of aluminum using a numerically controlled 
mill. This transition duct was also used as a mold to produce an identical fiber- 
glass duct used by Davis and Gessner for measurements made at the University of 
Washington [18]. 
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Table HI ! Transition duct parameters for Equation M.3 


x/D 

y o/D 

r,/D 

n y 

zq/D 

r z /D 

n z 

0.00 

0.00000 

0.500000 

2.00000 

0.00000 

0.5000000 

2.00000 

1.00 

0.00000 

0.500000 

2.00000 

0.00000 

0.5000000 

2’00000 

1.06 

0.00000 

0300165 

2.10467 

0.00000 

0.499854 

2.10467 

1.12 

0.00000 

0.501236 

2.21536 

0.00000 

0.498904 

2.21536 

1.18 

0.00000 

0.503911 

2.33308 

0.00000 

0.496531 

2.33308 

1.24 

0.00000 

0.508676 

2.45877 

0.00000 

0.492306 

2.45877 

1.30 

0.00000 

0.515822 

2.59342 

0.00000 

0.485967 

2.59342 

1.36 

0.00000 

0.525474 

2.73825 

0.00000 

0.477408 

2.73825 

1.42 

0.00000 

0.537602 

2.89471 

0.00000 

0:466651 

2.89471 

1.48 

0.00000 

0.552048 

3.06449 

0.00000 

0.453840 

3.06449 

1.54 

0.00000 

0.568540 

3.24972 

0.00000 

0.439214 

3.24972 

1.60 

0.00000 

0.586718 

3.45288 

0.00000 

0.423092 

3.45288 

1.66 

0.00000 

0.606151 

3.67704 

0.00000 

0.405858 

3.67704 

1.72 

0.00000 

0.626356 

3.92613 

0.00000 

0.387938 

3.92613 

1.78 

0.00000 

0.646823 

4.20500 

0.00000 

0.369787 

4.20500 

1.84 

0.00000 

0.667028 

4.51993 

0.00000 

0.351867 

4.51993 

1.90 

000000 

0.686461 

4.87910 

0.00000 

0.334633 

4.87910 

1.96 

0.00000 

0.704640 

5.29319 

0.00000 

0.318511 

5.29319 

2.02 

0.00000 

0.721131 

5.77751 

0.00000 

0.303885 

5.77751 

2.08 

0.00000 

0.735577 

6.35245 

0.00000 

0.291074 

6.35245 

2.14 

0.00000 

0.747705 

7.04838 

0.00000 

0.280317 

7.04838 

2.20 

0200000 

0.757357 

7.91061 

0.00000 

0.271758 

7.91061 

2.26 

0.00000 

0.764503 

9.01307 

0.00000 

0.265419 

9.01307 

2.32 

0.00000 

0.769268 

10.00000 

0.00000 

0.261194 

10.00000 

2.38 

0.00000 

0.771943 

10.00000 

0.00000 

0.258821 

10.00000 

2.44 

0.00000 

0.773015 

10.00000 

0.00000 

0.257871 

10.00000 

2.50 

0.00000 

0.773180 

10.00000 

0.00000 

0.257725 

10.00000 

5.47 

0.00000 

0.773180 

10.00000 

0.00000 

0.257725 

10.00000 
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Swirl Generator 

Rotating flow 

The function of the swirl generator was to create a solid body rotational flow 
superimposed on a uniform, one-dimensional flow. In cylindrical coordinates this 
velocity field is described by Equation m.5, where ^centerline * s ve I° c ity at the 
swirl generator centerline and 0 is the angular velocity of the solid body rotation. 


% = Centerline. C = flr, V. = 0 (HI.5) 


An example of this flow is represented in Figure m.5. The swirl angle, <f>, is given 
by Equation HI. 6. This is the angle of the flow relative to the x-axis measured by a 
stationary observer. This angle increases as the radius increases but it is V$, not <f), 
that increases linearly with the radius. 


<f> = arctan 



(HI.6) 


Ideally, this rotating flow would be free of wakes or other disturbances caused by 
the mechanism which creates the rotation. In principle this can be achieved by flow 
through a sufficiently long pipe rotating about its axis at a constant angular velocity. 
For the desired case of no radial velocity, the radial momentum equation requires 
that the radial pressure gradient balance the centrifugal forces acting on the fluid, as 
shown in Equation III.7. 


dp 


V? 


dr ^ r 


(DI.7) 


This condition is referred to as simple radial equilibrium [48]. For solid body rotation, 
where V$ — fir, the radial pressure gradient is given by Equation 113.8. Equation III. 8 
requires the static pressure distribution to increase parabolically in the radial direction 
from the centerline. For the rotating pipe, the Coriolis force is negligible since the 
axis of rotation corresponds to the principal direction of fluid motion. 
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Z 



Figure III.5 Solid body rotation superimposed on uniform flow 

| l = (HI. 8) 

Solid body rotation and the resulting static pressure distribution is centrifugally 
stable in the sense that a fluid particle whose position is perturbed radially will tend to 
return to its original position by the influence of the radial static pressure distribution. 
This stability is unlike, the case of flow between concentric cylinders when the inner 
pipe is Totaling and outer pipe is stationary [49] . 

, /fir 

<t> — arctan I — 

V *'* 

The amount of swirl for this flow is given by Equation III.9. The amount of swirl 
is determined by the incoming velocity and the angular velocity, which may be varied 
independently. This flow field will have a constant axial component of vorticity 
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of magnitude 2Q everywhere, and negligible radial and tangential components of 
vorticity everywhere except within the boundary layer. 

In important consideration in the analysis of a rotating pipe is whether work 
is added or extracted from the flow. Ignoring the boundary layer for the moment, 
consider an initially irrotational, incompressible, inviscid, uniform flow entering a 
rotating pipe. The flow entering the pipe will have constant radial total and static 
pressure distributions. The radial static pressure distribution of the exiting flow will 
increase parabolically from the centerline at a rate determined by the angular velocity 
of the rotating flow in order to maintain radial equilibrium. If no work is added or 
extracted from the flow the radial total pressure distribution of the exiting flow will 
remain constant and at the same level as the entering flow. This will require the total 
velocity magnitude and hence the axial component of velocity to decrease radially for 
the exiting flow. If work is added to the flow, like in a compressor, the exiting flow 
will have a radial total pressure distribution that will also increase parabolically from 
the centerline. The radial decrease of the axial component of velocity of the flow 
exiting the pipe will be less than the case when no work is added or extracted from 
the flow. If work is extracted from the flow, like a turbine, the radial total pressure 
distribution of the exiting flow will decrease parabolically from the centerline. The 
radial decrease of the axial component of velocity of the exiting flow will be greater 
than the case when no work is added or extracted from the flow. 

Swirl generator design 

A rotating pipe alone as a swirl generator was not satisfactory for this research 
for several reasons. A developing internal flow was desired. When shear stress has 
transferred momentum from the rotating pipe wall into the initially irrotational flow 
and solid body rotation is achieved the flow would be fully developed as well. If 
honeycomb were used at the entrance of the rotating pipe, like the swirl generator of 
Weske and Sturov [26], solid body rotation can be achieved before the flow is fully 
developed. However, the power requirements necessary to achieve solid body rotation 
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become excessive for large angular velocities and particularly for large through flow 
velocities. It was doubtful whether this energy could be transferred to the flow by 
the honeycomb without destroying the honeycomb. 

Another possible design was a set of stationary blades designed to create solid 
body rotation. This design would also create a rotating flow field immediately, but 
it was not satisfactory because of the presence of unwanted wakes created by the 
stationary blades and the difficulty in designing blades that exactly produce solid 
body rotation. 

Figure III.6 shows a schematic of the swirl generator that was conceived, designed 
and constructed for this study. This generator employs both stationary blades and a 
rotating pipe. The flow entering the stationary blades is nominally uniform and one- 
dimensional. The stationary blades create a near solid body rotation at an angular 
velocity that depends on both the axial flow velocity and the stationary blade design, 
which may be varied independently. Ideally, the angular velocity of the rotating pipe 
is given by Equation 131.10. 




V x tan <p 


(CLIO) 



Figure I II. 6 Swirl generator schematic 
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Immediately downstream of the stationary blades is the rotating pipe. The 
honeycomb and screen are fixed to the rotating pipe and these components rotate 
as one assembly. The honeycomb serves three functions: 

1. To serve as a rotor to impart motion to the rotating pipe. Ideally, the rotating 
pipe will achieve an angular velocity that matches the solid body rotation of the 
flow exiting the stationary blades. 

2. To reinforce the solid body rotation initiated by the stationary blades. The 
honeycomb is performing its customary roll by straightening the flow, this time' 
within the rotating frame of reference. 

3. To dissipate the wakes created by the stationary blades. 

A screen is located downwind of the honeycomb to dissipate the wakes created by 
the honeycomb. The final section of the swirl generator is the stationary exit section. 
This section functions as a stationary component to which the transition duct may 
be attached. 


Stationary blade design 

The terminology used to discuss blade design, shown in Figure III.7, has been 
adopted from Glassman [50]. Several assumptions were made to simplify the design 
of the stationary blades. It was presumed that the incoming flow is uniform and that 
along blade surfaces the flow satisfies the tangency condition. 

To accommodate the assumed incoming flow conditions the blade inlet angle, ipi, 
must be zero everywhere along the blade leading edge. The assumptions require the 
flow exiting the stationary blades be parallel to the camber line at the trailing edge. 
Therefore, the exiting flow is at an angle V>2 relative to the x-axis and the amount of 
swirl will be equal to the blade exit angle, <p = To achieve solid body rotation, 
the blade exit angle must satisfy Equation Ifl.ll. 


= arctan tan ^2, max) 


(fll.ll) 
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Tangent to camber line 



Figure UL7 Stationary blade design nomenclature 


In Equation III. 11 R is maximum radius and ^ 2 , max is the blade exit angle at r = R. 
The angular velocity of the flow is given by Equation HI. 12. 


n = 


V x tan ^2, max 

R 


(III. 12) 


This simple analysis cannot be used to deteraiine the shape of the blade from leading 
to trailing edge or the total number of blades required. 

A novel approach to blade design and construction was used to produce blades 
that exactly satisfy the blade inlet angle requirement and closely approximate the blade 
exit angle requirement, Equation HI. 11. These blades were fabricated by cutting the 
triangular shape a-b-e from a thin circular tube, as shown in Figure DI.8. The portion 
of the triangle used to construct the blade is the shaded shape a-b-c-d, the lower 
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Figure III.8 Blade fabrication from circular tubing 


part of the triangle, c-d-e, is truncated to allow attachment of the blade to a small 
annular center body. The blade leading edges were made blunt, the trailing edges 
were sharpened. 

Figure in.9 shows the blade design interpreted geometrically as the intersection 
of three circular cylinders and two planes. Cylinders 2 and 3 are concentric, with 
their centerlines along the z-axis. Cylinder 2 represents the inside of the inlet pipe, 
cylinder 3 represents the annular center body. Plane 2 is equivalent to the y 2 -plane, 
plane 1 is equivalent to an inclination of the yz-plane about the y-axis by the angle 
6 . The centerline of cylinder 1 lies in plane 1 and is normal to the y-axis. Cylinder 
1 contains the blade surface which is shaded, the dark shaded surface is visible, the 
light shaded surface is hidden. Each blade edge is the intersection of two surfaces, 
each blade comer the intersection of three surfaces. The blade edges, comers, and 
their corresponding surfaces are listed in Table III.2. 
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Plane 1 



Figure III.9 Geometric interpretation of blade design 
Parametric equations which describe the surface of the three cylinders are: 
Cylinder 1 : 

C\ : oi,o;2 — ► z,y,2 oil > 0, —7 r < a2 < 7r 
x = P sin 02 cos 8 — <*i sin 0, y = P — P cos o ;2 (HI. 13) 

z = P sin c*2 sin 8 a. \ cos 0 


Cylinder 2: 


C 2 : flu h — > x,y, 2 /9i > 0, -7T < fa <* 
x - —flu y = R sin/32, 2 = jRcos /?2 


(Til. 14) 
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Table III.2 Blade comer and edge construction from geometric surfaces 


Comer 

Surfaces 

Edge 

Surfaces 

a 

Cylinder 1 
Cylinder 2 
Plane 1 

a-b 

Cylinder 1 
Cylinder 2 

b 

Cylinder 1 
Cylinder 2 
Plane 2 

b-c 

(trailing edge) 

Cylinder 1 
Plane 2 

c 

Cylinder 1 
Cylinder 3 
Plane 2 

c-d 

Cylinder 1 
Cylinder 3 

d 

Cylinder 1 
Cylinder 3 
Plane 1 

d-a 

(leading edge) 

Cylinder 1 
Plane 1 


Cylinder 3: 


C?, : 71,-72 -» x, y, z 71 > 0, -7r < 72 < 7r 
x = -71, y = p sin 72, z = p cos 72 


(m.i5) 


To design and construct the blades required calculating the location of the 
blade comers. Each comer was determined by finding the intersection of the three 
corresponding surfaces. The Cartesian coordinates of point a is x = —R tan 6 , y = 
0, z = R. To find the location of point b. Equation HI. 16, a quadratic form in cos 02 , 
was solved to find 02 - This result was used in Equation III.17 to find a\ and the two 
are substituted into Equation in. 13 to determine the Cartesian coordinates. 


cos 2 a 2 + 2 tan 2 8 cos 8 + tan 2 6 



(HI. 16) 


Q] = P sin Q'2 cot 8 


(III.17) 
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Alternatively, Equation III. 18, a quadratic form in sin /? 2 , was solved to find /%. 
This result and = 0 were substituted into Equation HI. 14 to verify the Cartesian 
coordinates of point b. 


s “ 2 *- 2 (l) 


sec 2 8 -J- tan 2 6 = 0 


(m. 18 ) 


The Cartesian coordinates of point d are x = — ptantf, y = 0, z = p. The 
location of point c was found by both the procedures used to find point b, by either 
substituting p for R in Equation 111,16 or substituting p for R in Equation in. 18 and 
72 for 02 in IIL 18 and using III. 15 rather than HI. 14 to find the Cartesian coordinates 
of point b. 

With the Cartesian coordinates of a point on the blade surface known, the blade 
camber angle was determined from Equation HI. 19. 


ip 2 = arctan 


vV + z 2 ( x cos 9 + z sin 9) 
zP sec 9 + xy sin 9 — yz cos 9 


(HI. 19) 


Only two parameters detemiine the design of this type of stationary blade, P and 
9. The location of the four comers and blade exit angle were calculated for different 
parameter values and compared with two criteria in mind: 

1. The maximum blade exit angle, ip2, max> equal 20°. This value was consistent 
with swirl angles reported in studies of isolated single and multiple stage turbine 
flow exit angles [19-21]. 

2. The chord length of the blade be reasonably long. This is an attempt to minimize 
adverse pressure gradients along the blade surface. 

The values used for construction are P = 10.16 cm and 9 — 21.5°. The total 
number of blades used was 24. The consequence of too many blades (excessive wakes 
and total pressure loss) or too few blades (single airfoil behavior rather than cascade 
behavior) were considered, but the determination of the number of blades used was 
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very subjective. Table HI.3 lists the Cartesian and cylindrical coordinates of the blade 
constructed. Figure HI. 10 shows how the blade exit angle of the designed blade, given 
by Equation HI. 19, deviated from the ideal blade exit angle, given by Equation ELI 1. 
Table HI.3 Cartesian coordinates of stationary blade 
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CHAPTER IV. EXPERIMENTAL 
MEASUREMENTS AND INSTRUMENTS 

This chapter ^teseribes^the experimental measurements tmade,^ ^ 
techniques used, and the instruments employed during the course of this research. 
The Trace Gas Measurements section describes the trace gas technique in consider- 
able detail. The instruments and technique used for trace gas measurements were* 
developed specifically for this project and represent one of the major components 
of fins dissertation, in -addition, an improve^^^ calibra- 

tion and data reduction was developed. This procedure is referred to in the Pressure 
Measurements section, but the complete description has been deferred to Appendix 
A because of its length. 

Surface Oil Film Visualization 

A surface oil film visualization technique was used to obtain preliminary infor- 
mation about the flow field near the surface of the transition duct. Surface oil film 
visualization involves applying a mixture of oil and pigment to the duct surface. The 
air flowing over the surface induces motion within the oil film, depositing a streak of 
pigment indicating the direction of the flow field near the surface. This initial “quick 
look” was performed before other investigations and the results were used to guide 
the detailed aerodynamic -and trace gas measurements. 

Validity of surface oil film visualization 

Surface oil film visualization is best for observing steady flows due to the slow 
response of the oil film. Transient phenomena such as those associated with flow 
start up are not followed by the oil. Two important questions regarding the validity 
of surface oil film visualization, even for steady flows, are: 
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1. Does the presence of the oil film alter the boundary conditions of the flow to 
such a degree that it results in a different flow field? 

2. Does the pattern produced by the oil film accurately represent the flow field near 
the surface? 

An analysis by Squire etal. [51] showed the most significant parameter affecting 
the performance of surface oil film visualization is the ratio of the absolute viscosity 
of the boundary layer fluid to the viscosity of oil. For small viscosity ratios and thin 
oil films the boundary conditions are not materially changed by the presence of the oil 
film. Furthermore, for small viscosity ratios and large local shear stresses the direction 
indicated by the oil flow will accurately represent the flow field direction near the 
surface. Where the shear stress is small (e.g., near separation) the oil film motion 
becomes influenced by near surface pressure gradients and may not accurately reflect 
the near surface flow field direction. Squires analysis did not include gravitational 
forces acting on the oil film. In practice this can significantly alter the oil film pattern 
on other than horizontal surfaces in a manner similar to “runs” that appear in wet 
paint. Fortunately, these gravitational affects are relatively easy to identify. 

Surface oil film visualization technique 

The surface oil film visualization technique used is similar to methods described 
by Gessner and Chan [52] and Jurkovich et al. [53]. The resulting air to oil viscosity 
ratio was 10 -6 . A commercial dye which fluoresces when illuminated with ultraviolet 
light was mixed with the oil. Using fluorescent dye with ultraviolet illumination 
produces brighter images and reduces glare from the duct surface. 

At the beginning of a visualization test, the upper half of the transition duct 
was removed and the oil and dye mixture was applied to the lower duct surface. 
Care was taken in applying the oil and dye mixture to avoid a visible pattern. The 
upper half of the duct was replaced and the test facility was started, operated at test 
conditions for several minutes, then shut down. The upper half of the transition duct 
was subsequently removed and the lower half was illuminated with ultraviolet light 
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to excite the fluorescent dye. The resulting oil flow pattern was photographed with a 
35mm camera with sufficient exposure time to use ultraviolet light for illumination. 

Pressure Measurements 

Surface static pressure measurements 

Surface static pressure measurements were made through small (0.51 mm) tap 
holes in the transition duct surface. The tap hole axes were oriented normal to the 
duct surface. There were 50 surface static pressure taps equally spaced along the duct 
surface in the x z-plane. "Hie location of ffie smfacem^ 

in Figure IV. 1 as a solid line along the lower surface of the duct. The equipment 
used for measuring and recording pressure is described in the Manometry and Data 
Acquisition section of this chapter. 



Figure IV. 1 The location of surface static pressure measurements 
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Five-hole probe measurements 

A calibrated five-hole probe was used to measure flow field velocity direction, 
total pressure and static pressure at numerous points. From the total and static 
pressures the corresponding Mach number was calculated. This value of Mach 
number and the flow direction defined the Mach vector. The procedure used for 
five-hole probe calibration and data reduction is described in Appendix A. 

The five-hole probe Figure IV.2 shows a typical five-hole probe. Each of the 
five tubes are connected to a pressure transducer. All five-hole probe measurements 
were made using the yaw nulling technique. In this technique the five-hole probe is 
rotated around the axis of the probe shaft until the pressure measured at tubes two 
and four are equal. The angle through which the probe was rotated, the yaw nulling 
angle, was recorded along with the five pressure measurements. The equipment used 
for measuring and recording these pressures is described in the Manometry and Data 
Acquisition section of this chapter. 



Figure IV.2 A typical five-hole probe 
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Figure IV.3 The location of the four cross stream aerodynamic measurement planes 

Five-hole probe measurement locations Five-hole probe access to the flow 
was provided by machined openings in the transition duct. The axis of each opening 
was located to lie within one of the four cross stream aerodynamic measurement 
planes shown in Figure IV.3. For the designation of the measurement planes the 
letter suffix A in Figure IV.3 refers to aerodynamic and is used to differentiate the 
aerodynamic measurement planes from the trace gas measurement planes. The trace 
gas measurement planes, which are described in the next section, carry the letter suffix 
T and have axial locations that are different than for the aerodynamic measurement 
planes. 

For each aerodynamic measurement plane, the number of probe openings, their 
location, and the orientation of their axes varied. Figure IV.4 shows each probe 
opening axis in each measurement plane as a broken line. The first plane had nine 
probe openings located symmetrically about z-axis. The axis of each probe opening 
passed through the duct centerline. There were thirteen probe openings in the second 
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Figure IV .4 The location of the transition duct probe opening axis 
plane and fifteen probe openings each in the third and fourth planes. In these planes 
the probe openings were located symmetrically about the 2 -axis, and their axes were 
parallel to the 2 -axis. The middle probe opening of each plane passed through the 
duct centerline. 

Because the transition duct is symmetric with respect to the xy- and x 2 -planes 
it would have been sufficient to make measurements in only one quadrant of each 
measurement plane for flow without inlet swirl. Measurements were made in two 
quadrants, however, on both sides of the x 2 -plane for flow with and without inlet 
swirl. These additional measurements were made to confirm that the flow field does 
reproduce the symmetry of the transition duct for the case without inlet swirl, and to 
provide a comparison for the inlet swirl case, where flow field measurements in two 
adjacent quadrants are necessary. Approximately 480 five-hole probe measurements 
were recorded at each measurement plane for each test condition. 
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Five-hole probe positioning equipment The five-hole probe was positioned 
by a mechanism with five degrees of freedom. Three degrees of freedom were 
involved in positioning the probe at the probe opening and placing the probe shaft axis 
in alignment with the probe opening axis. This positioning was controlled manually 
.and remained fixed while measurements were made at that probe opening. The 
remaining two degrees of freedom were controlled by an electromechanical actuator 
and included translation along the probe shaft axis to position the probe for the 
measurement, and rotation of the probe about the probe shaft axis to null the probe.* 

Pitch and yaw angle reference measurements The procedure of making a set 
of measurements at one probe opening is referred to as a survey. The data reduction 
procedure used requires a reference measurement to establish pitch and yaw reference 
angles to which all angles measured during a survey are referred. One method of 
obtaining this reference is to make the reference measurement in a different facility 
with known flow field properties, usually a free jet. This procedure requires physically 
moving the actuator and probe from one facility to the other. In practice it is difficult 
or impossible to exactly reproduce the orientation of the actuator and probe between 
different facilities. 

Early trace gas tests in the transition duct suggested no deviation in the air flow 
along the duct centerline from the centerline. This result established the air flow along 
the duct centerline as a direction reference. The five-hole probe positioning equipment 
was designed to move the actuator and probe to any survey location without having to 
remove the actuator and disturb its orientation with respect to the positioning system 
or the transition duct. Before each survey a reference direction was established at the 
duct centerline. This procedure avoided the potential source of errors that exist when 
the reference direction measurements are made in a different facility. 

Survey procedures For each survey the actuator and probe were first moved to 
the survey probe opening. The actuator was translated to its fully extended position 
with the probe free of the actuator and held stationary to avoid contact with the 
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transition duct. Then, the probe was slid through the opening until it rested on the 
lower surface of the transition duct. The probe was then secured to the actuator, 
establishing the orientation between the probe and the actuator. This procedure 
established the reference for actuator position measurements and guaranteed that 
the probe could not be damaged by the actuator translating the probe into the duct 
surface. Once this orientation was established it was not disturbed until the survey 
was complete. The probe was then moved to the center probe opening of the plane to 
make the direction reference measurement. Then, the probe and actuator were moved 
back to the survey probe opening where the survey measurements were made. 

Preliminary measurements indicated the surface boundary layer was nominally 
7.62 mm thick. Each survey began at the surface of the lower half of the duct and 
proceeded upwards toward the horizontal plane of symmetry. Thirteen measurements, 
equally spaced at 0.63 mm intervals, were made in the first 7.62 mm of each survey. 
Measurements were continued at 2.54 mm intervals until the probe had crossed the 
zy-plane. All measurements were made in the duct half opposite the probe opening to 
avoid aerodynamic disturbances near the duct surface created by the probe opening. 

Trace Gas Measurements 

The objectives of the trace gas technique were described earlier, in Chapter II. 
This section focuses on the equipment and procedures used to make the trace gas mea- 
surements. The equipment necessary to make aerodynamic trace gas measurements 
can be divided into three major components: 

1. The trace gas injection system. 

2. The trace gas sampling system. 

3. The trace gas concentration measuring system. 

Because there are no commercially available trace gas equipment items specifically 
intended for aerodynamic tests like those presently described, the equipment must be 
assembled from components designed for other applications. 
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Trace gas injection system 

The trace gas injection system was designed to introduce the trace gas into the 
flow field at a point and with a specific flow rate as unobtrusively as possible. The 
injection system included an injection probe, a probe positioning system, a trace gas 
flow Tate control system, and a trace gas supply system. 

Trace gas injection probe Three outcomes considered important in the design 
of the injection probe are: 

1. Minimal downstream flow disturbance from the probe. 

2. An injection probe rugged enough iooperaie in the free stream at Mach numbers 
as great as 0.5. 

3. An injection probe capable of injecting trace gas at a sufficient rate to allow de- 
tection throughout the downstream flow field using the flame ionization detector. 

Some disturbance of the flow downstream of the injection probe is unavoidable. 
One obvious disturbance minimization step is to make the injection probe as thin 
as possible. However, this may lead to some pitfalls. The flow field downstream 
of the injection probe may be reasonably approximated by considering a circular 
cylinder in cross flow. Roshko’s [54] investigation of the wake of a circular cylinder 
in cross flow suggested that a stable, periodic, Karman vortex street occurs for a 
range of Reynolds number, based on cylinder diameter, from 40 to 150. This wake 
street decays very slowly and maintains its periodic structure for a long distance 
downstream. For a Reynolds number in the range from 300 to 10,000 the vortex 
street is quickly dissipated and by 50 diameters downstream the wake is completely 
turbulent. Therefore, the desire to make the probe as thin as possible was tempered 
by the requirement to keep the Reynolds number greater than 300, as well as less than 
10,000. For the velocity levels anticipated in the transition duct tests, for M = 0.50, 
this required the probe diameter to be greater than 0.03 mm and less than 1.04 mm. 

The review of trace gas literature in Chapter II briefly described trace gas injection 
methods. Several of the studies cited in Chapter II included the results of experiments 
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designed to measure and minimize disturbances from L-shaped type probes. Smith 
[29] and Moore [30] compared ethylene concentration distributions produced from 
trace gas injected through 3.2 mm and 1.1 mm outside diameter L-shaped probes into 
the boundary layer within a linear turbine cascade. TTiey concluded the trace gas was 
probably affected very little by probe diameter. The 3.2 mm probe was used for their 
turbine cascade measurements. 

Gallimore [32] and Gallimore and Cumpsty [33] compared three different injec- 
tion probe designs in a uniform free stream and within a turbulent boundary layer over 
a flat plate. The first probe was constructed from a straight length of 1.0 mm outside 
diameter tube, closed at one end with a 0.2 mm opening drilled through the tube wall 
near the closed end. The other two injection probes had an outside diameter of 0.56 
mm but differed in design, one was L-shaped and the other a “crooked” L-shape with 
an offset opening. The difference in concentration distributions resulting from trace 
gas injection by the different probes was most apparent in the free stream experiments. 
The smaller probes produced noticeably less trace gas spreading. The concentration 
distributions produced by the two smaller probes were nearly identical. The straight 
probe and the crooked L-shaped probes were used for the compressor measurements. 

Moore et al. [42] tested the performance of several different probes in a uniform 
free stream, a turbulent boundary layer, and a linear turbine cascade. Their free stream 
results showed the smaller probe, 1.07 mm, produced concentration distributions with 
less spreading, indicating the probe itself was contributing to turbulent mixing of the 
trace gas. In the boundary layer and linear turbine cascade, however, the effect of 
the probe size was much less apparent. Taken together, these two results indicate that 
in regions of low turbulence intensity some minor probe disturbance is unavoidable, 
while in regions of higher turbulence intensity, the disturbance created by the probe is 
not of significant magnitude compared to the existing turbulence to materially effect 
the mixing of the trace gas. 

All the trace gas injection probes described in the cited literature were operated 
in comparatively low speed flows and utilized a cantilever design. For the injection 
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probe to be as small as possible and yet withstand the aerodynamic loading imposed 
when it is extended to the center of the free stream at Mach numbers as great as 
0.50 required abandoning the cantilever design of previously used probes. Instead, a 
catenary design was used as shown in Figure IV.5. The injection tube is supported 
along its full length by a second tube which is held in tension across the rigid frame. 
The cross section of the rigid frame and the wall bounding the flow are shown in 
Figure IV.5. The catenary probe design requires an additional probe opening in the 
wall bounding the air flow. 

The trace gas injection probe was constructed from stainless steel hypodermic 
needle tubrng with an outsi^ -diameter of 0.56 mm. The catenary design allowed the 
trace gas injection probe to be as slender as the smallest reported in the literature, 
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yet be operated at significantly higher flow field velocities than is capable with a 
cantilever configuration. 

The flow rate of trace gas was chosen so that the bulk velocity of trace gas at the 
injection probe exit would approximately match the free stream velocity. The area of 
the trace gas injection probe exit was 0.067 mm 2 . For the transition duct experiments, 
with the free stream Mach number nominally equal to 0.50, this required a trace gas 
volumetric flow rate of approximately 600 ml/min. 

TVace gas injection probe positioning system Trace gas was injected from the 
seven different locations shown in Figure IV.6. All trace gas injection locations were 
within an axial plane situated one radius upstream of the beginning of the transition 
duct. One injection location was coincident with the centerline. The remaining six 
injection locations were placed radially at r/R = 0.498 and r/R = 0.995 along the 
^-axis and along lines through the centerline 60° on either side of the 2 -axis. 



Figure IV.6 Trace gas injection locations and measurement planes 
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Figure IV.7 Trace gas injection probe positioning equipment 


The equipment used to position the trace gas injection probe is shown in Figure 
IV.7. The trace gas injection probe was positioned by controlling two degrees of 
freedom, the radial movement and the circumferential movement. The radial position 
was controlled by moving radially the rigid frame which holds the trace gas injection 
probe in tension. A stepper motor and potentiometer were used to control and measure 
this radial motion. The rigid frame was attached to a large bearing whose center 
of rotation was coincident with the transition duct centerline. The circumferential 
motion was manually controlled by rotating the rigid frame around its rotation axis 
and accessing the flow through different pairs of probe access holes. 


Trace gas flow rate control system To achieve consistent trace gas results 
required carefully controlling the rate at which trace gas was injected. This was 
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accomplished in the transition duct test by using a Datametrics type 825 mass flow 
controller. This equipment allowed the desired mass flow rate of trace gas to be 
specified in engineering units at the controller panel, and it maintained the desired 
flow rate to within 1% of the specified flow rate. The controller measured the flow 
rate by heating a fraction of the gas and sensing the resulting temperature increase. 
Details of the principles and operation of the equipment are contained in Reference 55. 

The mass flow controller was calibrated by timing the displacement of water in an 
inverted graduated cylinder by gas passing through the controller. Since the pressure- 
and temperature of the gas were known, the mass of the gas was calculated from the 
volume of displaced water. This procedure was repeated several times for several 
different mass flow rates, and this information was used to adjust the controller. 
Because of the heat transfer principle used by the controller to determine the flow 
rate, the calibration results were, in general, only valid for a specific trace gas. 

Trace gas supply system. In this study, ethylene was used as the trace gas. 
Ethylene possesses several properties which make it desirable for trace gas studies: 

1. The molecular weight of ethylene, 28.1, is nearly the same as air, 29.0. This 
eliminates the effect of buoyancy forces on the motion of the trace gas. 

2. Under normal circumstances, ethylene and other hydrocarbon gases are only found 
in very minute quantities in air. Thus, the problem of having to adjust the trace 
gas concentration measurements to compensate for the background level of trace 
gas in air is avoided. 

3. The mean concentration of ethylene, or any hydrocarbon gas, can be measured 
very accurately with a flame ionization detector. This instrument is described 
later in this section. 

A potential hazard of ethylene is its flammability. The lower and upper flamma- 
bility limits of ethylene in air are 31,000 ppm and 320,000 ppm, respectively, and 
the autoignition temperature is 763 °K [56]. For the transition duct experiments, the 
injection rate of ethylene was such that far downstream, where the ethylene is well 
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mixed with air, the resulting concentration was approximately 2 ppm, well below 
the lower flammability limit. Undoubtedly, there are locations within the flow field 
where the concentration levels are within the flammability limits, but the flow field 
total temperature was nominally 294 °K, well below the autoignition temperature. 

The supply system included a supply cylinder of ethylene, a two-stage pressure 
regulator at the cylinder, and the associated tubing. Because of the flammable nature 
of ethylene, a small supply cylinder was used and the cylinder was stored outdoors 
when not in use. All ethylene supply connections were leak checked each day before 
the experiment. 

Trace gas sampling system 

The purpose of the trace gas sampling system was to provide the trace gas 
concentration measurement instrument with a steady and continuous supply of flow 
field sample from a specified location. The components of the trace gas sampling 
system included the sampling probe or tube, the probe positioning system, and the 
sample pump and tubing necessary to maintain the flow of sample from the sampling 
probe to the trace gas concentration measurement instrument. 

Trace gas sampling probe Trace gas sampling probes that have been used in 
previous studies were reviewed in Chapter II. For the references cited in Chapter II 
the measurements were performed at relatively low velocities, in nearly all cases the 
Mach number was less than 0.10. One difficulty in applying the trace gas technique 
to higher speed flows has been the .assumed need of isokinetic sampling of the trace 
gas [29], sampling at a velocity through the probe opening that matches the local host 
flow velocity. Isokinetic sampling is actually not necessary to accurately measure the 
local concentration of trace gas in a host gas flow. Injecting trace gas into the flow 
field creates an Eulerian concentration field. The local value of the concentration 
in the flow field does not depend on the rate at which the sample is acquired. The 
only possible pitfall resulting from the sampling rate is if the velocity through the 
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probe opening greatly exceeds the local host flow velocity, in which case the spatial 
resolution ability of the probe deteriorates. 

Another sampling issue is the importance of aligning the sampling probe precisely 
in the direction of the host flow. This concern is also largely unnecessary if the 
proper precaution is taken. Because of the measurement technique, the concentration 
of the sample will not depend on the orientation of the probe when the sample 
is acquired, provided the presence of the probe doesn’t materially alter the flow 
field. The sampling probe is usually the same size as a pressure probe, where the 
disturbance caused by the probe is usually acceptable. These principals of sampling 
rate and sampling probe alignment were confirmed during the initial tests of the trace 
gas instruments and technique described in Chapter V. 

If a conventional trace gas sampling probe design had been employed in the 
present study, probe access would have been through the same openings used for 
five-hole probe access that were described in the Pressure Measurements section of 
this chapter. However, the distance between these openings, 1.91 cm, was too large 
to satisfactorily resolve the concentration gradients expected in the yz- plane. 

A novel method for sampling was designed to allow finer spatial resolution in the 
yz-plane while still using the existing five-hole probe access openings. The design 
employed as the sampling probe a flexible tube that extended upstream to the sampling 
plane while firmly base mounted at a fixed location at the exit of the transition duct. 
The tube was constructed from stainless steel hypodermic needle tubing with outside 
and inside diameters of 1.83 mm and 1.37 mm, respectively. The position of the tube 
opening was controlled by the trace gas sampling probe positioning equipment, and 
the flow field sample was drawn through the tube opening. 

TVace gas sampling probe positioning equipment The trace gas sampling 
probe and positioning mechanism are shown in Figure IV.8. The location where 
the probe is mounted at the duct exit lies downstream beyond the right border of 
this photograph. The position of the probe opening can be visualized as lying on 
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the surface of an imaginary scribed sphere whose radius is equal to the probe length. 
Because the probe length was sufficiently large, the deviation of the surface of the 
scribed sphere from the measurement plane was negligible. 

The position of the sampling probe arm was controlled by the same electro- 
mechanical actuator used to position five-hole probes. Both the translation, or axial, 
and rotation, or radial, modes of the actuator were used to position the trace gas sam- 
pling probe opening to the desired position in the measurement plane. The location 
of the trace gas measurement planes were shown in Figure IV.6. Different length 
sampling tubes were used for each measurement plane. 

A computer program was written to relate the actuator radial and axial coordinates 
to transition duct coordinates. Figure IV.9 shows an idealization of the trace gas 
sampling probe and positioning system. The transformation is given by Equations 



Figure IV. 8 Trace gas sampling probe positioning equipment 
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IV. 1 . These coordinate transformations are an idealization, they represent a close 
approximation but fail to account for deviations such as deflections of the sampling 
probe tube. 

R(a,2 + r sin 6 — 62) 

y = 02 + 
z = 


z 



Figure IV.9 Parameters used for the trace gas 
sampling system coordinate transformation 
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To produce a closer approximation, a\, 02, 61,62, R and r were considered vari- 
ables rather than constants and the set of functions given by Equations IV. 1 were 
calculated while varying these six parameters. The error of Equations IV. 1 was de- 
fined to be the sum of the magnitude of the difference of the measured and calculated 
sampling tube opening positions. A numerical scheme known as the steepest decent 
method was used to search the parameter space to minimize this error. In using this 
method a beginning point was entered and the gradient of the error is approximated 
using finite differences. The negative of the gradient is a vector in the parameter 
space which points in the direction of values which result in the most rapid decrease 
of error. The search proceeded along this direction until the error no longer decreased 
at which time the gradient was again calculated. The search was then repeated in the 
new direction and this process was continued until the search in a new direction no 
longer produced any further decrease in the error. 

The typical procedure involved measuring the position of the sampling probe tube 
in transition duct coordinates for many different actuator axial and radial positions. 
The six physical quantities which correspond to the parameters were also measured 
and used as the starting condition of the search. The numerical method was then 
used to calculate the parameters which were used in the coordinate transformation 
program. There is no guarantee that this method resulted in parameter values that 
produced a global minimum in error, or even a local minimum. However, in all 
cases it did significantly reduce the error, the average error usually being only about 
0.25 mm. This procedure was repeated for each actuator location for each of the 
four measurement planes. 

Trace gas sampling pump and tubing A sampling pump was required to sup- 
ply the trace gas concentration measurement instrument with a continuous supply 
of sample. This pump was necessary to overcome the pressure loss in the tubing 
between the sampling probe opening and the concentration measurement instrument 
and because the total pressure of the flow field was less than atmospheric pressure. 
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An important consideration was whether the sampling pump should be located 
upstream or downstream of the gas concentration measurement instrument. Smith [29] 
and Gallimore [32] used a sampling pump upstream of the trace gas concentration 
measurement instrument, the sampling pump of Bauer [34] was located downstream. 
The advantage of a downstream location is the pump cannot contaminate the sample 
with hydrocarbons, which are present in the pump lubricating oil. The disadvantage 
of a downstream location is that much of the tubing from the trace gas sampling 
probe to the pump, including the trace gas concentration measurement instrument, 
will operate below atmospheric pressure. This allows the possibility of atmospheric 
air leaking into the sampling system and diluting the sample. When the sampling 
pump is placed upstream of the trace gas concentration measurement instrument the 
pressure downstream of the pump will be greater than atmospheric. Any leaks of 
the sample fluid out will not change the concentration of the sample. For this study, 
the sample pump was located upstream of the trace gas concentration measurement 
instrument. A K.N.F. Neuberger, Inc., Model N05ANI pump was used. This is an 
oil free pump designed specifically for contamination free pumping of gasses. 

Stainless steel tubing was used throughout the trace gas sampling system. A 
special cleaning procedure was used to remove any hydrocarbons from the tubing 
before installation. First, the tubing was continuously flushed with acetone for several 
minutes. This was accomplished with a pump and large reservoir of acetone. The 
pump forced the acetone through the tubing which discharged back into the reservoir. 
Then, this procedure was repeated with methyl alcohol to remove any acetone residue. 
Finally, the tubing was dried with a continuous flow of nitrogen for several minutes. 

Trace gas concentration measurement system 

The instrument used for is largely determined by the choice of trace gas. Re- 
searchers that use helium frequently use a thermal conductivity cell [27, 39, 40, 37]. 
Carbon dioxide concentration is commonly measured with an nondispersive infrared 
detector [57, 38, 37]. A flame ionization detector provides a very accurate measure- 



ment of ethylene concentration and has been used extensively for trace gas studies 
[28-35, 43, 58, 42], and was used for this investigation. 

Trace gas concentration detection equipment used in previous studies was briefly 
reviewed in Chapter II. A GOW-MAC Instrument Company Model 23—500 Total 
Hydrocarbon Analyzer was used for this research project This instrument utilizes 
a flame ionization detector for measuring hydrocarbon concentration. The flame 
ionization detector may be thought of as a carbon atom counter. It consists of a 
small burner within which fuel and oxidizer are mixed and burned. Hydrogen and 
zero-gas air (air containing less than 0.5 ppm total hydrocarbons), were used as fuel 
and oxidizer, Tespectively, so there are no hydrocarbons present m this flame. When 
the sample is mixed with the flame, any hydrocarbons are burned and the carbon 
atoms become ionized. The carbon ions and electrons pass between two electrodes, 
decreasing the resistance between the electrodes and thus permitting an electric current 
to pass. This current is directly proportional to the amount of carbon ions present. 
Additional details about the principles and operation of this instrument are provided 
in Reference 59. The electronic output of the total hydrocarbon analyzer was routed 
to the data acquisition system, described in the next section. 

The flame ionization detector responds to both the amount of carbon atoms present 
in the sample, and the rate at which the sample passes through the detector. This 
makes it very important to maintain a constant flow rate of sample through the flame 
ionization detector. This was accomplished within the total hydrocarbon analyzer 
with a “tee” connector in the sample supply line with a capillary tube restrictor on 
the leg of the tee to the flame ionization detector and a back pressure regulator on 
the other leg. This arrangement provides a constant flow rate of sample to the flame 
ionization detector regardless of the sample flow rate into the tee, provided, of course, 
the flow rate into the tee exceeded the desired flow rate to the detector. The fraction 
of sample fluid not directed to the flame ionization detector was discharged to the 
atmosphere. Therefore, fluctuations in the flow rate of the trace gas sampling pump 
did not affect the flow rate into the flame ionization detector. The suggested flow 
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rate of sample fluid into the flame ionization detector specified by the manufacturer, 
25 ml/min, was used. 

It is also important to maintain a constant flow rate of fuel and oxidizer to the 
flame ionization detector. Pressure regulators followed by capillary tube restrictors 
upstream of the flame ionization detector were used to maintain constant flow rates. 
The suggested flow rates of fuel and oxidizer specified by the manufacturer, 25 ml/min 
and 245 ml/min, respectively/ were used. Because of the extreme flammability of 
hydrogen, a small hydrogen supply cylinder was used. This cylinder was stored- 
outdoors when not in use, and all connections were leak checked daily during tests. 

The flame ionization detector was calibrated by preparing a mixture of ethylene 
and zero-gas air of known concentration and supplying this as a sample to the total 
hydrocarbon analyzer. This calibration was performed periodically throughout the 
period of trace gas testing. After the initial calibration, it was seldom necessary to 
make any further adjustments. 


Manometry and Data Acquisition 

Pressure measurements 

Five-hole probe pressures and surface static pressures were measured with a Pres- 
sure Systems Incorporated Model 780B electronically scanned pressure measurement 
system. The system is described in Reference 60. 

The system consists of four major components; the systems controller, the data 
acquisition and control unit, the pressure calibration unit, and the sensor modules. A 
schematic diagram of these components is shown in Figure IV. 10. Each experimental 
pressure was measured with an individual electronic pressure transducer contained 
in a sensor module. For these experiments, four sensor modules, each containing 
32 individual differential pressure transducers, were used. Each individual pressure 
transducer had a full scale range of 34.5 KPa, and was referenced to atmospheric 
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pressure so the difference between experimental pressures and atmospheric pressure 
was measured. 



Figure TV.10 Schematic of the Pressure Systems Incorporated 
Model 780B electronically scanned pressure measurement system 


The pressure calibration unit included pneumatic valves to switch the input of 
the individual pressure transducers from the experiment pressure to a calibration 
pressure. The calibration pressure was also measured with a single, high-accuracy, 
pressure transducer contained in the pressure calibration unit. For these experiments 
a Digiquartz absolute pressure transducer with a full scale range of 103.4 KPa and 
an accuracy of 0.02% of full scale, or 21 Pa, was used. 

The measurement of the calibration pressure by the Digiquartz transducer was 
used to calibrate each individual pressure transducer. The calibration procedure 
was controlled by the data acquisition and control unit. Three different calibration 
pressures, spanning the range of pressure anticipated in the experiment, were used to 
calculate separate calibration curves for each individual pressure transducer. These 
calculations were performed by the data acquisition and control unit. By calibrating 
the differential transducers in the sensor modules with an absolute transducer the 
resulting calibration allowed the differential transducers to measure absolute pressure. 
This calibration procedure was performed at the beginning of each test day and was 
repeated at approximately 30 minute intervals during tests. 
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Measurement of the individual pressure transducers and the frequency of their 
calibration were controlled by the data acquisition and control unit which received 
its commands from the systems controller, a personal computer. While acquiring 
data, the data acquisition and control unit measured each experimental pressure 
sequentially at a rate of 0.00001 seconds per measurement. This was repeated 
until eight measurements were recorded for each transducer. The arithmetic mean of 
these measurements was calculated and converted to an absolute pressure using the 
calibration. This pressure was then passed to the systems controller. 

Data acquisition 

The purpose of the data acquisition was to display and record all of the experi- 
mental data as well as ancillary data which describe the facility operating condition. 
This task was performed by an existing data acquisition hardware and software pack- 
age called Escort 2, which was developed at NASA Lewis Research Center. The 
data acquisition system accepted both digital input, such as the pressures measured 
by the electronically scanned pressure measurement system, and analog input, such 
as the voltage output of the total hydrocarbon analyzer or the potentiometers used 
to measure probe position. 

The data were continuously updated, converted to engineering units, and displayed 
on a CRT monitor in the control room of the test facility. When requested, the data 
were recorded by a centrally located minicomputer. Archival files of all recorded 
data were maintained. The data files were transferred from the minicomputer to a 
computer engineering workstation for data reduction, display and analysis. 

Hot Wire Measurements 

Thermal anemometry 

A schematic of the hot wire measurement equipment is shown in Figure IV. 1 1 . 
All measurements were made with a single normal wire probe. A standard TSI model 
1260A-T1.5 hot wire probe with a 4 /im diameter tungsten sensing element was used. 
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A TSI model IFA100 constant temperature anemometer was used to acquire hot wire 
measurements. The details of this anemometer are described in Reference 61. The 
anemometer bridge voltage was output to a Hewlett Packard model 7603 oscilloscope, 
a Fluke model 8920 RMS voltmeter and a Masscomp model MC 5450 computer and 
high speed digital data acquisition system. All hot wire measurements were recorded 
with die digital data acquisition system. The oscilloscope and RMS voltmeter were 
used only to detect problems at run time. The analysis of the digitized hot wire data 
were performed with Sun computer workstations. 

Attempts to make hot wire measurements using a standard TSI hot wire probe 
holder were unsuccessful for Mach numbers greaterthan 0. 10 because of aerodynamic 
buffeting of the probe holder. This buffeting was detected by the oscilloscope 
representation of the signal which became less random and more periodic as the 
Mach number increased. A much stiffer probe holder with approximately the same 



Figure IV. 11 Schematic of hot wire anemometry equipment 
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cross stream projected area as the TSI probe holder but having a highly elliptical 
cross section was fabricated in order to make hot wire measurements at the higher 
velocities. Both the oscilloscope and frequency analysis of the digitized signal showed 
this new probe holder worked satisfactorily at velocities up to the highest Mach 
number recorded, 0.50. 

Hot wire calibration 

To calibrate the hot wire probe, the probe was placed at the centerline of the 
straight pipe in the Internal Fluid Mechanics Facility with the wire oriented normal 
to the mean flow direction. The mean flow velocity was varied and determined with 
total and static pressure measurements. The mean anemometer bridge voltage and 
the mean velocity were recorded and used for calibration. Data were recorded at 
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Figure IV. 12 Hot wire calibration results 
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five overheat ratios, 1.2, 1,4, 1.6, 1.8, 2.0. The calibration data were approximated 
with a King’s law [62] relation. Equation IV.2. A least squares procedure was used 
to determine the value of the constants B and n for each overheat ratio. The value 
of A was determined from zero velocity measurements. Figure IV. 12 shows the hot 
wire calibration results. The data recorded at die lowest overheat ratio, 1.2, were 
unsatisfactory and is not presented. 

E 2 = A + B(pU) n (IV.2) 
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CHAPTER V. RESULTS 

The data presented in this chapter can be divided into two categories; data gath- 
ered to validate an experimental facility or experimental technique, and data acquired 
to examine the aerodynamics of the transition duct flows. The first three sections 
present validation data for the trace gas technique, the Internal Fluid Mechanics Fa- 
cility and the swirl generator. The last two sections present the results of the transition 
duct measurements for flow without and with inlet swirl. Analysis of the transition, 
duct results is contained in Chapter VI. 

All of the results presented in this chapter are nondimensional. Total and static 
pressures are presented as pressure coefficients Pq and p*, having been nondimen- 
sionalized as indicated in Equations V.l and V.2. Measured velocity vectors were 
divided by the local sonic velocity c to yield a Mach vector. This result was then 
normalized by the centerline Mach number, as shown in Equation V.3, resulting in 
a velocity coefficient M*. Bold type in Equation V.3 indicates vector quantities. In 
Equations V.l through V.3 the reference variables, subscripted centerline (r = 0) or 
wall {r = R), were evaluated at a location one radius upstream of the transition duct 
inlet ix/D = -0.5). This location represents conditions near the exit of the straight 
pipe or swirl generator or the inlet of the transition duct. 


P o 


Po ~ Pwall 
Po, centerline — Pwall 


(V.l) 


P Pwall 

Po, centerline — Pwall 


(V.2) 


M* = 


V/c 


(V.3) 


-^centerline 

The concentration measured in trace gas surveys was the ethylene mole fraction, 
given by Equation V.4. In Equation V.4 Ethylene re f ers to the molar amount of 
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ethylene present in a sample and n refers to the total molar amount of the quantity 
of matter in the sample. The results are presented as an ethylene concentration 
coefficient C*, where the local ethylene concentration has been normalized by the 
peak value of ethylene concentration measured for that particular distribution, as 
shown in Equation VS. 


^ "ethylene 

^ethylene — ~ 


(V.4) 


c * _ ^ethylene 
^ethylene,peak 

Plotting contours of the trace gas and aerodynamic data presented a challenge 
because the locations of the data points were not on a regular grid. The irregular 
grids occurred when there was an unequal number of data points in the surveys in 
a measurement plane or when additional data were taken in a region to resolve a 
flow feature. Similar problems resulting from irregularly distributed data occurred 
in the statistical analysis of the trace gas results. No available software could 
satisfactorily Tesolve both these problemssoanumerical computer program was 
written to interpolate the irregularly distributed data onto a regular grid. A portion 
of this program was also used for statistical analysis of trace gas data. The details of 
the interpolation procedure are presented in Appendix B. 


Trace Gas Method Validation 

Prior to conducting the transition duct experiments, preliminary trace gas mea- 
surements were made in the calibration tunnel described in Chapter in to insure 
all the trace gas equipment was operating properly. These measurements were also 
performed to reveal information about the trace gas technique itself. To test the 
spatial resolution of the sampling probe, ethylene was injected at the centerline of 
the calibration tunnel jet and the flow was extensively sampled in a cross stream 
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plane approximately 60 mm downstream of the injection point. Figure V.l shows a 
contour plot of the measured distribution of the concentration coefficient, C*. This 
figure shows that the trace gas system is capable of resolving flow features as small 
as approximately 1.0 mm. 

A test was performed to determine the sensitivity of the ethylene concentration 
measurements to the ethylene injection rate. Ethylene was injected at the center- 
line of the calibration tunnel jet Two separate cross stream sample surveys were 
conducted approximately 110 mm downstream of the injection point. The rate of 
ethylene injection varied between the two surveys while all other test conditions were 
unchanged. A comparison of the measured concentration distributions, shown in Fig- 
ure V.2, when presented as concentration coefficients, show very little difference. 
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Figure V.l Concentration coefficient C* distribution for preliminary 
tests to ascertain the spatial resolution of the trace gas technique 
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Figure V.2 Concentration coefficient C* distributions for 
preliminary trace gas tests to ascertain the sensitivity of 
concentration measurements to the ethylene injection flow rate 

The rate of ethylene injection for the result shown on the right hand side of Figure 
V.2 was 1.5 times greater than that used for the left hand side result. Within some 
reasonable limitations, the concentration coefficient appears to be insensitive to the 
rate at which ethylene is injected into the host flow. 

To test the sensitivity of the ethylene concentration measurements to a mismatch 
between the sample velocity at the probe opening and the iiost flow velocity, two cross 
stream sample surveys were compared. For both surveys ethylene was injected at the 
centerline of the calibration tunnel jet and the flow was sampled in cross stream planes 
approximately 110 mm downstream of the point of injection. For the results shown 
on the left hand side of Figure V.3 the jet Mach number was approximately 0.10. 
The jet Mach number was approximately 0.60 for the results shown on the right hand 
side Figure V.3. The rate at which the sample was drawn was held constant, insuring 
a significant velocity mismatch between these two cases. The results presented in 
Figure V.3, however, show little difference in the concentration coefficient. These 
results demonstrate that isokinetic sampling is not necessary. 

To test the sensitivity of the ethylene concentration measurements to the sampling 
probe alignment, an L-shaped type sampling probe was employed. This probe was 
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Figure V.3 Concentration coefficient C* distributions for preliminary trace gas tests 
to ascertain the sensitivity of concentration measurements to the sampling flow rate 

constructed so the probe opening was in line with the probe shaft axis so that the 
probe could be rotated, changing the probe orientation, without changing the location 
of the probe opening. This allowed the probe to be intentionally misaligned with 
the host flow direction without changing the probe opening location. Ethylene was 
injected at the jet centerline and measurements were made downstream at a single 
location. The probe alignment was displaced approximately 30° either side of the host 
flow direction without any discernible change in the measured sample concentration. 
This evidence suggests that at least within this range of angles, the concentration 
measurements are insensitive to misalignment of the probe with the host flow. 

Internal Fluid Mechanics Facility Validation 

Measurements were made of the velocity distribution and free stream turbulence 
near the straight pipe exit at a location one radius upstream of the transition duct 
inlet (x/D = —0.5). The purposes of these measurements were to assess the flow 
quality' of the Internal Fluid Mechanics Facility and to provide inlet conditions for 
CFD calculations of the transition duct flow field. 




y mm 
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Velocity measurements 

The velocity measurements were made within the straight pipe along a radial 
line from the wall to the centerline using a three-hole pressure probe. Data were 
recorded at approximately 0.64 mm intervals within the boundary layer and at 2.54 
mm intervals from the boundary layer edge to the centerline. Measurements were 
made at three different centerline Mach numbers, 0.20, 0.35, and 0.50, corresponding 
to the three nonswirling inlet conditions at which transition duct measurements were 
recorded. 

Figure V.4 shows the measured distribution of the axial component of the velocity 
coefficient M* at x/D = —0.5. Distances were measured from the centerline and 
normalized with the pipe radius. Outside the boundary layer the axial component 
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Figure V.4 The radial distributions of the axial component 
of the velocity coefficient M* at x/D = —0.5 without swirl 
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of velocity was notably uniform. There were no measurable tangential velocity 
components. Radial velocity components were assumed to be negligibly small. 

For each of the three centerline Mach numbers the boundary layer thickness 60.95 
(Equation V. 6 ), displacement thickness 8* (Equation V.7), momentum thickness 8 
(Equation V. 8 ), and shape factor H (Equation V.9) were determined by numerical 
integration of the experimental data. These results are summarized in Table V.l. 
Equations V.7 and V .8 have been modified from their conventional formulations by a 
factor of r/R within the integrand because of the axisymmetric geometry of the pipe.- 
Table V.l shows that the boundary layers are quite thin; in all cases the boundary 
layer thickness is less than 10 % of the pipe radius and the displacement thickness is 
less than 2% of the pipe radius. The shape factors measured indicate the boundary 
layers are turbulent and there is little streamwise pressure gradient. 

^ = sup jl — — \U(r) < 0.95l7 center ]j ne | (V.6) 
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~R 



U \ / r \ / r \ 

^centerline/ 


(V.7) 
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^centerline 
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^centerline 



(V.8) 


8 * 

H=j (V.9) 

Figure V.5 shows the velocity data within the boundary layer plotted in law of 
the wall coordinates. The conventional formulations of the law of the wall coordinate 
variables (Equation V.10) were used. The distance y refers to the distance from the 
wall, y = R — r. The friction velocity, u*, defined by Equation V.ll was determined 
by calculating the value of u* that provides the best approximation of the experimental 
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Table V.l Boundary layer thickness #o.95> displacement thickness 8*, 
momentum thickness 9 and shape factor H at x/D = —0.5 without swirl. 


-^centerline 

£o. 95/-R (Xl00) 

8*/R (xlOO) 

9/R (xlOO) 

H 

0.20 

6.85 

1.37 

0.95 

1.42 

0.35 

8.01 

1.49 

1.03 

1.44 

0.50 

7.52 

1.41 

0.99 

1.42 


data by Equation V.12, where k = 0.41 and B = 5.5. These results also indicate no . 
noticeable anomalies in the boundary layer entering the transition duct. 
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Figure V.5 The axial velocity distribution in law of 
the wall coordinates at x/D = —0.5 without swirl 
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(V.ll) 

u+ = j]ny + + B (V.12) 

lYirbulence measurements 

Hot wire measurements were made at the centerline of the straight pipe section 
of the Internal Fluid Mechanics Facility for flow without swirl. Measurements were 
made for Mach numbers ranging from 0.05 to 0.50 and at five different overheat ratios. 
These data were gathered to ascertain the turbulence intensity Tu and turbulence 
integral length scale A in the free stream flow. The intention of these measurements 
was not to study the turbulence in detail but rather to determine velocity and length 
scale values which in general terms characterize the turbulence. 

The calibration of the hot wire probe, described in Chapter IV, and the acquisition 
of the turbulence data were conducted simultaneously. The known mean flow 
information was used to determine the calibration parameters while the digitized signal 
was analyzed to discover the turbulence information. A 50 Khz data acquisition 
sampling rate was used and the anemometer output was recorded for 2 seconds, 
acquiring at total of 100,000 recorded anemometer bridge voltage values for each 
Mach number and overheat ratio. The mean and RMS values of the bridge voltage 
were calculated numerically. The turbulence intensity was calculated using Equation 
V.13 for each Mach number and overheat ratio. These results indicate turbulence the 
intensity is nominally 0.65%. 




fl(r) = V(i )"l i ±- r ) 


(V.14) 
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The definitions of the turbulence integral scales and their significance are de- 
scribed by Bradshaw [63]. To determine the integral length scale required first calcu- 
lating the autocorrelation R(t ) of the velocity fluctuations given by Equation V.14. 
Rough estimates indicated the value of the integral time scale was no larger than 
approximately 0.015 seconds. Since data were recorded for 2 seconds the whole 
data signal was partitioned and 20 independent values of the autocorrelation were 
calculated. An ensemble average of all 20 autocorrelations was then calculated. This 
average autocorrelation was then integrated to determine the integral time scale T 
as shown in Equation V.15. Taylor’s hypothesis was then invoked to use the mean 
velocity to calculate the integral length scale by Equation V.16. Ibis process was re- 
peated for each Mach number and overheat ratio. The results showed some variation, 
but the nominal value of the integral length scale was approximately 3 mm. 

OO 

T = J R(r)dT (V.15) 

0 

a = Ft (v.i6) 

Swirl Generator Validation 

Measurements were made of the operating conditions (e.g., centerline Mach 
number and maximum ideal swirl angle), velocity and pressure distributions near the 
exit of the swirl generator at a location one radius upstream of the transition duct inlet 
(x/D = —0.5). This data provided information about the operating characteristics of 
the swirl generator as well as inlet conditions for CFD calculations of the transition 
duct flow field. 

Swirl generator operating condition 

To determine the swirl generator operating conditions and for nondimensional- 
izing experimental results, knowing the centerline total and static pressures near the 



71 


swirl generator exit and the angular velocity of this flow was essential. A series 
of experiments was therefore conducted to explore how these operating conditions 
could be determined. 


Centerline conditions During transition duct experiments a way to determine 
the centerline total and static pressure at x/D = —0.5 without using an intrusive 
probe was sought in order to avoid flow field disturbances in the transition duct. 
Without the swirl generator the centerline total pressure was assumed to be equal to 

the upstream settling chamber pressure (Po, centerline ~ ^settling chamber) the 
centerline static pressure was taken to be equal to the static pressure measured at the 

pipe wall (Pcenterline == Pwall) - The wa ^ sta ^ c pressure was measured with two 
surface static pressure taps located at x/D = —0.5. The centerline Mach number 
was calculated using Equation V.17. 


M. 


centerline — 




Po, centerline \ 


\ W Pcenterline / 


- 1 


7-1 


(V.17) 


Two critical assumptions in this procedure are no total pressure loss between the 
settling chamber and inlet centerline and no radial static pressure gradient at the inlet. 
These are reasonable assumptions for flow without the swirl generator. However, 
neither assumption is valid with the swirl generator installed. There are significant 
total pressure losses associated with the blades, honeycomb, and screen of the swirl 
generator and a radial pressure gradient exists downstream of the swirl generator to 
balance centrifugal forces established in the fluid by the swirling flow. 

A series of experiments were conducted with the swirl generator installed to 
determine the relationship between the centerline total and static pressures at x/D — 
-0.5 and the settling chamber and wall static pressures. A calibrated five-hole 
pressure probe was used to measure the centerline total and static pressures while 
the swirl generator operating conditions were varied. The relationships between 
centerline total and static pressures at x/D = —0.5 and the settling chamber and 
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Figure V.6 The relationships between centerline total and static pressures 
at x/D = —0.5 and the settling chamber and wall static pressures with swirl 


wall static pressures are shown in Figure V.6. Different symbols mark experimental 
measurements recorded on different days. The solid line is a plot Equation V.18, the 
broken line is a plot of Equation V.19. 
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(V.19) 
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Both constants of Equations V.18 and V.19 were determined by a least squares 
linear approximation. Equations V.18 and V.19 allowed the centerline total and static 
pressures at x/D = —0.5 to be determined from the settling chamber and wall static 
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pressures, which were measured without an intrusive probe. The ratio of total to 
static pressure at the centerline was used in Equation V.17 to determine the centerline 
Mach number. 

Angular velocity A nonintrusive proximity sensor was used to measure the 
angular velocity of the swirl generator. This proximity sensor, which detects ferrous 
metal, was attached to the stationary outer casing. A steel screw was attached to 
the aluminum rotating pipe of the swirl generator. A voltage spike was produced^ 
in the output signal of the proximity sensor with each revolution of the aluminum 
pipe. A frequency counter was used to measure the time between consecutive 
voltage spikes. Figure V.7 shows the relationship between the swirl generator 
angular velocity and centerline Mach number. Different symbols denote experimental 



Figure V.7 The relationship of swirl generator 
angular velocity to centerline Mach number 




74 


measurements recorded on different days. The solid line shows a linear approximation 
of the data for 0 < < 0.4. For > 0.4 the data deviates 

significantly from the linear approximation. The standard deviation of the data 
for 0 < M centerline < 0.4 is 8.6 rad/s. At M centerline = 0.35, the operating 
condition used for the transition duct measurements with swirl, this equates to an 
angular velocity uncertainty of 2.8%. The maximum ideal swirl angle calculated 
using Equation III.9 is <f>max. = 15.6°. 

Swirl generator aerodynamic survey 

Measurements were made with a calibrated five-hole probe of the velocity, total 
pressure, and static pressure distribution near the exit of the swirl generator at 
x/D = —0.5. These data were used to assess the flow field produced by the swirl 
generator and provide inlet conditions for CFD calculations of the transition duct 
flow field. 

Velocity measurements Figure V .8 shows the distributions of the axial, radial 
and tangential components of the velocity coefficient M* at x/D = —0.5. The results 
are for three different centerline Mach numbers. 

The Tadial component of the velocity coefficient confirms that there is no mea- 
surable radial velocity component at the swirl generator exit. The tangential velocity 
coefficient component data show that tangential velocity increases nearly linearly with 
radial distance from the centerline, achieving the desired objective of solid body rota- 
tion. The axial component of the velocity coefficient has a small aberration near the 
centerline. This is a result of the total pressure losses produced by the small annular 
center body used to hold the stationary swirl generator blades. 

The boundary layer thickness 60 . 95 , displacement thickness 6 *, momentum thick- 
ness 6 and shape factor H were calculated from this data. The results are contained 
in Table V.2. The boundary layer entering the transition duct is thicker for the inlet 
swiTl case compared to the nonswirling case. The additional boundary layer thick- 
ness results primarily from work extracted from the flow to rotate the swirl generator 
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Figure V.8 The radial distributions of the axial, radial and tangential 
components of the velocity coefficient M* at xJD = —0.5 with swirl 

pipe against the bearing friction. The shape factors in table V.2 are higher than the 
nonswirling case. This is also probably a feature of the axial velocity deficit resulting 
from work extraction rather than an indication of an axial pressure gradient. 

Figure V.9 shows the velocity data plotted in law of the wall coordinates. Because 
the boundary layer is thicker than the nonswirling case the values of y + for the data 


Table V.2 Boundary layer thickness £o.95» displacement thickness 6*, 
momentum thickness 9 and shape factor H at x/D = —0.5 with swirl. 


-^centerline 

£ 0 . 95 /^ (xlOO) 

8*/R (xlOO) 

S/R (xlOO) 

H 

0.19 

19.99 

4.02 

2.41 

1.67 

0.26 

20.18 

4.03 

2.40 

1.68 

0.36 

20.48 


2.37 

1.71 
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in Figure V.9 are lower than in Figure V.5. There is some deviation from the law of 
the wall behavior at the lowest values of y + measured but, in general, the agreement 
appears quite good. 



Figure V.9 The axial velocity distribution in law 
of the wall coordinates at x/D = —0.5 with swirl 


Pressure measurements Figure V.10 shows the radial distributions of the total 
and static pressure coefficients measured at x/D — —0.5. The parabolic static 
pressure profile predicted by Equation III. 8 is apparent in Figure V.10. The losses in 
total pressure associated with the swirl generator annular center body near r/R = 0.15 
and from work extracted to rotate the swirl generator pipe against the bearing friction 
at r/R > 0.75 is conspicuous. 
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Figure V.10 The radial distributions of the total pressure coefficient 
Pq and static pressure coefficient p* at x/D = —0.5 with swirl 

The Transition Duct Without Inlet Swirl 

Table V.3 is a summary of the test conditions and measurement methods em- 
ployed for the transition duct experiments without inlet swirl. The aerodynamic mea- 
surements include surface static pressure measurements and detailed surveys with a 
five-hole probe in four cross stream measurement planes. For trace gas measure- 
ments, ethylene was injected at seven different locations in the injection plane and 
the trace gas concentration levels were measured in four cross stream measurement 
planes. The nondimensionalized aerodynamic measurement data vary little between 
test conditions, showing independence of inlet centerline Mach number or Reynolds 
number within the range of each considered. 
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Table V.3 Summary of test conditions and experimental 
methods for transition duct studies without inlet swirl 


■^centerline 

Rep 

Measurements recorded 

0.20 

918,000 

Aerodynamic 

0.35 

1,547,000 

Aerodynamic 

0.50 

2,086,000 

Aerodynamic, trace gas 



and surface oil film 



visualization 


Surface oil film visualization 

The first experiments conducted were surface oil film visualization in the transition 
duct without inlet swirl. Figures V.ll and V.12 show the flow visualization results 
from two different views. In both figures the flow is from the left to the right. A 
local coordinate system is shown in each figure to help establish the orientation of 
the photographs. These photographs were digitized with an image scanner and then 
numerically enhanced to improve their contrast. 

Trace gas measurements 

In Figures V.13 through V.16 the cross section of the duct is drawn to the 
same scale. Each figure shows contours of the ethylene concentration coefficient 
C* measured in that plane for all seven injection locations. The contour levels shown 
in Figures V.13 through V.16 are at 20, 40, 60, and 80 percent of the peak value. 
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Figure V.ll Surface oil film visualization without inlet swirl 
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Aerodynamic measurements 

Surface static pressure measurements Open symbols in Figure V.17 represent 
values of the static pressure coefficient p* distributed along the lower surface of the 
duct for all three test conditions. Shown as vertical broken lines are the locations of 
the four aerodynamic measurement planes. The solid symbols on these lines represent 
the static pressure coefficient measured at the duct centerline with a five-hole probe. 
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Figure V.17 Surface static pressure coefficient p* for flow without inlet swirl 

Five-hole probe measurements The results of the detailed five-hole probe 
surveys in four cross stream measurement planes within the transition duct are shown 
in Figures V.18 through V.41. The cross section of the duct is drawn to the same 
scale in each figure. Figures V.18 through V.25 show results at -^centerline = 0-0, 


• t i i i i t i 'i i | i ii i < i i i i ' j i i rr i . r i i i | i i " i i t r i i i . ' f i- v t t i i " i~i < | i > » " r ~i i r ~i i 



Plane 1A Plane 2A Plane 3A Plane 4A 

I I I 1 .1. 1 l JL L 1 1 t .1, 1 1 1...J i — 1 — 1 I I ..1. -L Ki I I I I t I I t 1 I I, 1— L-l— I .lilt i — UJ 1.1 I I I K .1 1 1 J i 





Figure V.18 Axial (upper) and transverse (lower) components of the velocity 
coefficient M* in plane 1A for M center jj ne = 0.20 without inlet swirl 




Figure V.19 Axial (upper) and transverse (lower) components of the velocity 
coefficient M* in plane 2A for ^centerline = 0-20 without inlet swirl 







Figure V.21 Axial (upper) and transverse (lower) components of the velocity 
coefficient M* in plane 4A for M center jj ne = 0.20 without inlet swirl 
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Figure V.26 Axial (upper) and transverse (lower) components of the velocity 
coefficient M* in plane 1A for ^centerline = 0-35 without inlet swirl 



Figure V.27 Axial (upper) and transverse (lower) components of the velocity 
coefficient M* in plane 2A for jlf C enterline = 0.35 without inlet swirl 








coefficient 




:oeffi< 






Figure V.34 Axial (upper) and transverse (lower) components of the velocity 
coefficient M* in plane 1A for M center ii ne = 0.50 without inlet swirl 



Figure V.35 Axial (upper) and transverse (lower) components of the velocity 
coefficient M* in plane 2A for -Renter line = 0-50 without inlet swirl 
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Figures V.26 through V.33 show results at = 0.35, and Figures V.34 

through V.41 show results at -^centerline = 0.50. 

Figures V.18 through V.21, V.26 through V.29, and V.34 through V.37 show 
all three components of the velocity coefficient M*. Vector plots of the transverse 
velocity coefficient components are shown in the lower half of the duct In the upper 
half of the duct are contour plots of the magnitude of the axial velocity coefficient 
component Beneath the duct in each figure is a vector labeled -^centerline* This * s 
the reference length scale used for the vector plots of the transverse velocity coefficient 
components. The same reference length was used in all figures. This allows direct 
comparisons ofthe velocity coefficientcomponemstobemadebetweenmeasurement 
planes. 

Distributions of the total pressure coefficient pq and the static pressure coefficient 
p* are shown in Figures V.22 through V.25, V.30 through V.33, and V.38 through 
V.41. Drawn in the lower half of the duct are contour plots of total pressure 
coefficient. Contour plots of the static pressure coefficient are drawn in the upper 
half of the duct. The same contour levels are used in all cases. The dashed lines that 
appear in some static pressure coefficient contours indicate negative values. 

In all cases the actual measurements were made in the lower half of the duct 
only, the data shown in the duct upper half are a rotation from the duct lower half. 

The Transition Duct With Inlet Swirl 

Surface oil film visualization and aerodynamic measurements were acquired 
for flow in the transition duct with inlet swirl. The aerodynamic measurements 
consist of surface static pressures and detailed aerodynamic surveys in four cross 
stream measurement planes. All transition duct measurements for the inlet swirl 
case were made at one condition, -^centerline = 0-35. The nonswirling flow 
aerodynamic data showed insignificant sensitivity to Mach number level within the 
range 0.20 < -^centerline ^ 0-50 and Reynolds number within the range 918,000 < 
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Re < 2,086,000. The locations of the surface static pressure measurements and 
the aerodynamic measurement planes are indentical to those for the nonswirling 
aerodynamic measurements. Trace gas measurements were not repeated for the 
swirling flow case. 

Surface oil film visualization 

Figures V.43 and V.44 show the results of the surface oil film visualization tests 
for the transition duct flow with inlet swirl. The same procedures used to digitize 
and enhance the nonswirling surface oil film visualization photographs were applied 
to these results. A local coordinate system is shown to help establish the orientation 
of these photographs. 

Aerodynamic measurements 

Surface static pressure measurements In Figure V.42 the transition duct sur- 
face static pressure coefficients p* distributed along the lower surface of the duct with 
inlet swirl are plotted with open symbols. The surface static pressure coefficients for 
flow without inlet swirl at the same centerline Mach number are also plotted for 
comparison. Shown as vertical broken lines are the locations of the four aerody- 
namic measurement planes. The solid symbols represent the centerline static pressure 
coefficients measured during the aerodynamic surveys. 

Five-hole probe measurements The results of the aerodynamic surveys of the 
transition duct flow with inlet swirl are shown in Figures V.45 through V.52. The 
same conventions used in plotting the nonswirling results have been used in the figures 
for swirling flow. The velocity plots are shown in Figures V.45 through V.48, the 
total and static pressure plots are shown in Figures V.49 through V.52. 




Figure V.42 Surface static pressure coefficient p* for flow with and without inlet swirl 
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Figure V.43 Surface oil film visualization with inlet swirl 
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Figure V.44 Surface oil film visualization with inlet swirl 
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CHAPTER VL ANALYSIS 

This chapter discusses the results of the measurements of the transition duct flow 
field. The focus of the first section is an analysis of the trace gas results. The 
statistics of the trace gas concentration distributions are presented and discussed. 
The procedure used to calculate the trace gas statistics is described in Appendix 
C. Aerodynamic measurements of the transition duct flow without inlet swirl are 
used to help interpret the trace gas results. The second section concentrates on the 
aerodynamic measurements of the flow through the transition duct and a comparison 
of the flow field with and without inlet swirl. In this chapter the nondimension al 
quantities Pq,p*, M* and C* are referred to simply as total pressure, static pressure, 
velocity and concentration. 


Trace Gas Results 

In a turbulent flow that is steady in the mean, the injection of trace gas will 
create a concentration field which itself is steady in the mean. Equation VI. 1 governs 
the transport of trace gas in a turbulent incompressible flow. The left hand side of 
Equation VI. 1 represents the convection of the mean concentration of trace gas by 
the mean flow, the first term on the right hand represents the molecular diffusion of 
the mean trace gas concentration, the second term represents the turbulent diffusion 
of the trace gas concentration. 


-dC —dC —dC o— 

U j— + v— + W — = D m V 2 C - 
ox oy oz 


( d^d) djv^d) 
y dx dy 


d(w'd) \ 
+ dz ) 


(VI. 1) 


For engineering calculations, mixing length theory is commonly used to relate the 
turbulent diffusion of any quantity (e.g., mass, momentum, energy) to the gradient of 
the mean value of the quantity, as shown in Equation VI.2. This creates a turbulent 
diffusion term in Equation VI. 1 that is similar to the molecular diffusion term. The 
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turbulent diffusivity coefficients are the product off a turbulent velocity and a turbulent 
mixing length, as shown in Equation VI.3. Often, the turbulent velocities are related 
to correlations of velocity fluctuations, which persist in the free stream, rather than 
gradients of the mean velocities, which vanish in the free stream. The mixing length 
is ameasure of ffie ctistance a fluid ^artcle ^ with the mean 

flow. In Equation VI.2, Dij is a turbulence velocity and length scale correlation 
coefficient defined by Equation VI.3. This creates a second order tensor in the 
transport equation for the rate of turbulent diffusion, shown in Equation VI.4, written 
in index summation notation [64]. 



When the trace gas technique is used to infer information about the flow field 
there is an ambiguity problem, referred to earlier in the literature review in Chapter 
II in the discussion of the paper by Wisler et al. [35]. The ambiguity problem arises 
when experimental measurements, which represent a solution to Equation VI.4 with 
a known boundary condition, are used to determine the coefficients of Equation VI.4, 
Ui and Dij, which are unknown. Although the solutions to Equation VI.4, with some 
restrictions on the coefficients and boundary conditions, are unique, Equation VI.4 
with different coefficients and identical boundary conditions can give rise to identical 
solutions. The ambiguity problem is exacerbated by the fact that the experimental 
data are not known everywhere, but only at a finite number of locations where 
measurements are made. This makes it difficult to interpret trace gas results without 
either making some assumptions which simplify Equation VI.4, or supplementing the 
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trace gas data with additional information, or both. The supplemental information 
may be experimental data or numerical results. 

For homogeneous turbulence, the turbulent diffusion tensor can be written in terms 
of a symmetric tensor. With the additional simplification of assuming the turbulence 
to be isotropic, the turbulent diffusion tensor reduces to a single turbulent diffusion 
coefficient Dt. Even when the turbulence is not homogenous and not isotropic, it is 
common in engineering calculations to use a single turbulent diffusion coefficient. 

The mixing length theory of turbulence is analogous to the kinetic theory of' 
gases where the molecular diffusivity coefficient is the product of velocity and length 
scales that govern molecular transport, the sonic velocity c, and the molecular mean 
free path A. For molecular diffusion, the diffusivity coefficient is essentially constant 
throughout the flow field, however, for turbulent diffusion, the coefficient changes 
with location. 


D t v t l t 

rv 

Djn c\ 


(VI.5) 


The relative importance of turbulent and molecular diffusion is expressed as the 
ratio of the two diffusivity coefficients, which varies with location. An estimate 
of the magnitude of this ratio is given by Equation VI.5. In the free stream, the 
turbulent velocity and integral length scale were determined experimentally, and the 
ratio is approximately 180. By analogy with the turbulent viscosity coefficient, in the 
boundary layer this ratio increases across the wake region and then decreases across 
the overlap region, approaching zero at the wall [65]. This is adequate justification to 
neglect the molecular diffusion of trace gas in further analysis. Neglecting molecular 
diffusion and using a single turbulent diffusion parameter reduces Equation VI.4 to 
Equation VI.6. 


+v& c +w ac = 

ox oy oz 


(VI.6) 
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Another assumption commonly used is to neglect diffusion in the stream wise 
direction. For the simple case of a uniform flow in the x direction, this changes 
Equation VI.6 from an elliptic to a parabolic partial differential equation, like the 
heat equation, where the x dimension corresponds to time in the heat equation. 

The following analysis will consider one speciaicaseoflhe^^ 
equation, when the molecular diffusion is negligible, the turbulent diffusion is ho- 
mogeneous but not necessarily isotropic, the turbulent diffusion in the axial direction 
is negligible, and the mean flow is uniform and in the x direction. Under these re- 
strictions, Equation VI.4 reduces to Equation VI.7. Without any loss in generality, 
it is assumed the y- and 2 -axis are coincident with the principal axis of the tensor 
Dij, so that D yz = 0. 


^d£ 

dx 


d 2 C 

,yy dy 2 


d fdC\ 

U ^~ = D «-nr + 2D V*fy (- ) + D ** dz 2 


d 2 C 


(VI.7) 


C = 


su 


47 T X . D yy D z 


6Xp— 4r 


U(y 2 z * 


+ 


yy 


D z 


(VI.8) 


2 Dyy2x 
cr — — = — 

y U 


D zz 2x 

U 


(V 1.9) 


For the case of a point source of trace gas located at the origin, the solution to 
Equation VI.7 is given by Equation VI.8. The variances of the solution are given 
by Equation VI.9, the covariance is zero because of the choice of axes. Although 
equation VI.7 appears simplified beyond the point of relevance, its solution provides 
some useful insight into interpreting the trace gas results. For this particular case: 


1. The contours of the solution will always be concentric ellipses 

2. The variances of the solution will always increase with axial distance 

3. The area of the trace gas spreading can be estimated by Equation VI. 10 which 
is the area of the ellipse defined by Equation VI. 11 
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(VI. 11) 


For each of the 28 concentration distributions measured, the mean values, vari- 
ances, and covariances of the concentration distributions were calculated. These 
values are presented in Tables VI. 1 through VIA Also presented in these tables are 
the covariance correlation coefficient and an estimate of the spreading area. The sta- 
tistics are useful to help quantify characteristics of the trace gas measurements, like 
the amount of trace gas spreading or the shape and orientation of the contours. The 
definition of these statistics and the numerical method used to calculate their values 
are described in Appendix C. All values in the tables are nondimensional. 

The locations of the four cross stream measurements planes differed between the 
trace gas and aerodynamic measurements, as was described in Chapter IV. This was 
a result of the difference in hardware required for each method. The locations of the 
measurement planes for each method are summarized in Table VI.5. 

It is useful to use the aerodynamic results to assist in interpreting the trace gas 
data when possible. The location of trace gas plane 3T and aerodynamic plane 2A are 
identical. Figure VI.2 shows the transverse velocity components from aerodynamic 
plane 2A superimposed over the results from trace gas plane 3T. Trace gas plane 
2T and aerodynamic plane 1A are near enough that it is reasonable to superimpose 
the transverse velocity components from aerodynamic plane 1A and contours from 
trace gas plane 2T. This is shown in Figure VI. 1. Trace gas plane 4T lies between 
aerodynamic planes 3A and 4A. These planes are located in the constant cross section 
region of the duct. Gradients in the axial direction are small enough that it is 
reasonable to interpolate aerodynamic measurements between planes 3A and 4A to 
compare with the trace gas measurements at plane 4T. This result is shown in Figure 
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Table VI.2 Plane 2T trace gas concentration distribution statistics 


129 




Table VI.3 Plane 3T trace gas concentration distribution statistics 
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Table VI.4 Plane 4T trace gas concentration distribution statistics 
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Table VI.5 Location oftrace gas andaerodynamic measurementplanes 


Trace gas measurements 

Aerodynamic measurements 

Plane 

X 

Plane 

X ‘ 


D 


D 

IT 


1A 

1.49 

2T 

1.43 

2A 

1.99 

3T 

1.99 

3A 

2.55 

4T 

3.36 

4A 

3.93 


VI.3. No comparison of trace gas measurements with aerodynamic measurements is 
possible at trace gas plane IT. 

For trace gas injected at the transition duct centerline, y/D = 0.0 z/D = 0.0, the 
covariance remained relatively small along the length of the duct This indicates that 
the turbulence field is symmetric along the duct centerline. These trace gas contours 
became increasingly elongated in the y direction at downstream locations. This can 
be seen in the ratio of Cy to a\. Two possible explanations for this are either the 
turbulence is anisotropic, that is y/v® > Vufi or the concentration distribution is 
being affected by the transverse components trf 'velocity. Two experimental facts 
suggest that the elongation is produced by the transverse velocity components. First, 
the value of a\ is decreasing from plane IT to plane 3T and then increases from plane 
3T to plane 4T. It would be impossible for o\ to decrease by turbulent diffusion alone 
in a uniform flow field. Second, velocity measurements show flow converging in the 
xz-plane and diverging in the zy-plane. The region of decreasing <r| corresponds to 
the region where transverse velocity components are greatest. This can be seen in 
Figures VI. 1 and VI.2. At plane 4T, Figure VI.3, the transverse velocity components 
near the centerline are negligible and has increased between planes 3T and 4T. 

For trace gas injected at y/D = 0.000 z/D = 0.249 (12 o’clock, halfway between 
the duct centerline and the wall) the results are statistically similar to the results for 
trace gas injected at y/D = 0.000 z/D = 0.000. However, at planes 3T and 4T the 




IVI centerline 

Figure VI. 1 Ethylene concentration distributions with transverse velocity components at plane 2T. 
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contour shapes appear to be noticeably different. The contour shapes are no longer 
elliptic, rather they are bean shaped, having only one plane of symmetry. The simple 
model of homogeneous turbulence in a uniform mean flow allows only elliptic shaped 
contours, which would indicate the contours result from turbulence inhomogeneity, 
or from nonuniformity in the mean flow. The mean flow certainly is nontmiform. 
However, there is nothing to indicate the turbulence is not homogeneous. Thus, the 
bean shape of the contours probably results from transverse velocity components. 
Also, the amount of spreading is approximately the same as for trace gas injected at 
y/D = 0.000 z/D = 0.000 which indicates roughly the same turbulence conditions. 

For trace gas injected at y/D — 0.000 z/D = 0.498 (12 o’clock, near the duct 
wall) the statistics of the measurements indicate a very large elongation in the y 
direction. The results of Klebanoff [65] show that the magnitude of the transverse 
turbulent velocity fluctuations, and turbulent viscosity coefficient, increase greatly 
in a turbulent boundary layer. In a turbulent boundary layer the turbulence is very 
inhomogeneous. The turbulence decreases rapidly away from the duct surface, in 
this case the 2 direction, while retaining its magnitude in directions parallel to the 
wall, specifically the y direction. This should lead to preferential diffusion in the 
y direction. However, turbulent diffusion alone is not exclusively responsible for 
the observed contours because a 2 z again decreases from plane IT to plane 3T. The 
total amount of spreading is roughly twice that observed of trace gas injected at 
y/D — 0.000 z/D = 0.000 and y/D = 0.000 z/D = 0.249. It appears that both 
turbulent diffusion and transverse velocities are responsible for the large transport of 
trace gas in the y direction displayed by the contours. 

For trace gas injected at y/D = -0.215 z/D — 0.124 and y/D = 0.215 z/D = 
0.124 (about 10 and 2 o’clock, between the duct centerline and wall) the amount of 
spreading at each of the four planes is comparable to the spreading of trace gas injected 
at y/D — 0.000 z/D — 0.000 and y/D = 0.000 z/D = 0.249. This indicates that 
in the free stream the turbulence is nominally homogeneous. Each contour becomes 
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more elongated as it travels downstream. Unlike the results for trace gas injected at 
y/D = 0.000 z/D = 0.000 and y/D = 0.000 z/D = 0.249, the contours change 
their orientation. This can be observed in the covariance values or the correlation 
coefficient. The major axis of the contours are no longer parallel to the y-axis, but 
are inclined to the y-axis. This probably results from the transverse velocities. 

Trace gas injected at y/D = —0.430 z/D — 0.249 and y/D — 0.430 z/D = 
0.249 (about 10 o’clock and 2 o’clock, near the duct wall) produced the most 
interesting results. At planes IT and 2T the statistics are nearly identical to the' 
statistics for trace gas injected at y/D = 0.000 z/D = 0.498, except for the 
orientation. At plane 2T, however, the contours appear visibly different, assuming a 
cusped teardrop shape. At plane 3T the statistics differ significantly from the statistics 
of trace gas injected at y/D = 0.000 z/D - 0.498. The spreading was approximately 
twice as great. This probably resulted from the transfer of boundary layer fluid, with 
its greater turbulence intensity and diffusivity, over a large area near the side walls by 
the axial vortices. At plane 4T the much greater spreading is obvious. The unusual 
shape results from vortices near the transition duct side wall, which are discussed in 
the following section. 


Aerodynamic Results 

Flow without inlet swirl 

The distribution of static pressures in the transition duct flow was generally 
attributed to the response of the flow field outside the boundary layer to the changing 
duct geometry. The change in cross section of the duct requires the flow to converge in 
the x£-plane and diverge in the xy-plane. The duct wall deflected the incoming flow, 
initially parallel to the x-axis, away from the x-axis in the xy-plane and towards the 
x-axis in the xz-plane. For the duct flow without inlet swirl, the curved streamlines 
initially generate a saddle-shaped static pressure distribution in the j/2-plane, with 
minimum static pressures near either end of the y-axis and maximum static pressures 
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near either end of the z-axis. This static pressure distribution can be seen at plane 1A, 
shown in Figures V.19, V.27, and V.35. The maximum pressures are also observable 
as the maximum value that appear in the surface static pressure plot shown in Figure 
V.14. 

Another saddle-shaped static pressure distribution in the yz-plane was developed 
further downstream by curved streamlines when the flow was forced by the duct 
wall back to a direction that is nominally parallel to the re -axis in the constant 
cross section region following the transition region. This static pressure distribution 
involved maximum static pressures near either end of the y-axis and minimum static 
pressures near eitta^end^of die z-axis- Cvidenee of this dismbutio in 

the static pressure distributions at plane 3 A, shown in Figures V.21, V.29, and V.37 
and as the minimum value of the surface static pressure plot shown in Figure V.14. 
Miau et al. [16, 17] and Davis and Gessner [18] also observed this reversal of the 
static pressure distribution in their surface static pressure measurements. 

The location of the maximum and minimum static pressures were reversed 
between planes 1A and 3 A by the change in streamline curvature between these 
planes. In order to affect this reversal of static pressure, there must be an intermediate 
location where the static pressure distribution in the yz-plane is nominally flat. The 
results show the intermediate location must be between planes 1A, at x/D — 1.49, and 
2A, at x/D = 1.99. The static pressure distribution at plane 2A, shown in Figures 
V.20, V.28 and V.36, already shows the reversal of the location of the maximum 
and minimum static pressures, although the gradients are less than those observed at 
plane 3A. The surface static pressure measurements made along axial cross sections 
in an identical duct by Davis and Gessner [18] show that the reversal occurs between 
x/D = 1.60 and x/D = 1.90. The overall level of static pressure was higher at planes 
1 A and 2A because of the increase in cross-sectional area at these measurement planes. 

The static pressure distribution at plane 4A, shown in Figures V.22, V.30 and 
V.38, is nominally flat. The slight gradients observed in the static pressure distribution 
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* 




probably resulted from the upstream influence of the discharge plenum at the transition 
duct exit, which was an abrupt enlargement. 

Because the static pressure distribution was established by the flow field outside 
the boundary layer where the momentum is greatest, the resulting static pressure 
distribution can produce significant turning of the flow within the boundary layer 
where the momentum is less. This is the source of skew-induced secondary flow 
referred to by Bradshaw [66]. This can be seen in Figures V.16, V.24 and V.32, where 
the greatest transverse velocity components appear near the duct surface towards the 
side walls. 


The explanation of skew-induced secondary flows in the proceeding paragraph 
has its foundation in the vorticity transport equation. The x component of the vorticity 
transport equation, neglecting viscous and Reynolds stress terms, is given by Equation 
VI. 12. The first term in the right hand side of Equation VI. 12 is referred to as the 


vortex stretching term, the next two terms are referred to as the vortex tilting terms. 


+ ^y~p, h w 2' 


(VI. 12) 


Consider aside from the transition duct flow the idealized case of a simple boundary 


layer in the xz-plane, U(y),V = 0, W = 0, with a cross stream pressure gradient, 
|f = Jg = 0, jf ^ 0, Equation VI.12 reduces to VI.13. The right hand side 
of Equation VI.13 arises from the second term of the right hand side of Equation 


VI.12. 


_ 1 dpdU 

dx pU dz dy 


(VI.13) 


The outward directed flow in the boundary layer near the side walls seen in the 
transverse velocity plots at planes 1A and 2A was the genesis of the two pair of 
counter-rotating side wall vortices that were clearly present at planes 3A and 4A. 
Although the static pressure distribution reversed orientation between planes 1A and 
2A, the second static pressure distribution failed to completely reverse the secondary 
fluid motion initiated by the first static pressure distribution. This is so because the 



140 


outward directed flow in the boundary layer near the side walls followed the duct 
surface around the corner and was redirected inwards along the side walls before 
the static pressure distribution had reversed its orientation. This can be seen in the 
surface oil film visualization results. Figure V.9. Along the side wall the flow is now 
approaching the xy-plane from top and bottom. Continuity forces the converging flow 
outward away from the side walls and the vortex pattern was established. The reversal 
in the orientation of the static pressure distribution only drove the vortices further 
away from the side walls. This can be seen by comparing the plots of transverse 
velocity components at planes 3A and 4A. 

Davis and Gessner [18] observed nearly identical behavior in the transverse 
velocity in their measurement planes downstream of the transition region in an 
identical transition duct. Patrick and McCormick [14, 15] observed similar side wall 
vortices at the exit plane of their aspect ratio six duct. The side wall vortices were not 
apparent at the exit plane of their aspect ratio three duct, although outward directed 
flow in the boundary layer near the side walls, which gives rise to axial vorticity, 
was observed at the exit plane. Miau et al. [16, 17] observed axial vorticity near 
the exit plane of their aspect ratio two transition duct which is opposite in sign to the 
dominant side wall vortices seen in Figures V.18, V.26 and V.34. 

Figures VI.4 through VI.7 show contours of axial vorticity in the upper half 
of the transition duct cross section for flow without inlet swirl. The axial vorticity 
was calculated from a finite difference approximation using the transverse velocity 
components which are shown in the lower half of the cross section. Positive vorticity 
(counter clockwise) is represented by solid lines, -negative vorticity (clockwise) is 
represented by broken lines. 

In Figure VI.4, at plane 1A, the initiation of the side wall vortices can be seen 
in the generation of positive axial vorticity in the upper left quadrant and negative 
axial vorticity in the upper right quadrant. This axial vorticity resulted from pressure 
gradient driven, skew-induced secondary flow in the boundary layer, as explained 
earlier. The sign of the axial vorticity shown in Figure VI.4 is correctly predicted by 




Figure VI.4 Axial vorticity at plane 1A for ^ ce nterline = 0.35 without inlet swirl 



Figure Vl.5 Axial vorticity at plane 2A for M center ii ne = 0.35 without inlet swirl 




Figure VI.6 Axial vorticity at plane 3 A for ^centerline = 0-35 without inlet swirl 
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Equation VI. 13. In Figure VI.5, at plane 2A, convection had carried the axial vorticity 
outwards towards the side walls. At plane 3A, shown in Figure VI.6, convection had 
carried the axial vorticity observed at plane 1A to the side walls and out of the 
measurement region. The new static pressure distribution, however, had generated 
axial vorticity which is opposite in sign from the axial vorticity observed at plane 1 A. 
This provides additional evidence that the second static pressure distribution doesn’t 
simply eliminate the axial vorticity produced by the first static pressure distribution 
because convection had carried the vorticity observed at plane 1A away from the 
region where the opposite signed axial vorticity is produced. In Figure VI.7 the axial 
vorticity observed at plane 1A is now readily apparent near the side walls, while the 
much weaker opposite signed axial vorticity has remained in nearly the same location 
as observed at plane 3A. 

Although there is no aerodynamic data very near the duct side wall, the surface 
oil film visualization provides evidence of a two additional pair of vortices in this 
region. This can be seen in surface oil film visualization results presented in Figure 
V.9. In the constant cross section region downstream of the transition region a bright 
line appears on the duct side wall approximately 0.15 radius from the xy-plane. 
Nearby lines exhibit asymptotic behavior, approaching the bright line from either 
side. Figure VI.8 is a sketch representing the pattern of secondary flow inferred from 
the surface oil film visualization. Numerical calculations at NASA Lewis Research 
Center using Navier-Stokes codes support this interpretation of the surface oil film 
visualization [67]. 

An effect of the vortices observed downstream of the transition region was the 
convection of fluid from the boundary layer into the free stream region. This produced 
distortion in the axial velocity and total pressure contours. This is particularly apparent 
at plane 4A, in Figures V.19, V.27 and V.35, where the high total pressure loss fluid, 
normally associated with the boundary layer, extends outwards from the side wall 
towards the duct centerline. 
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Figure VI-8 Sketch of secondary flow patterns inferred from surface 
oil film visualization for transition duct flow without Met swirl 


There are no results which indicate any region of separated flow in the transition 
duct. Davis and Gessner [18] made this same observation in their study of an identical 
transition duct. Miau et al. [16, 17] observed a region of separated flow along the 
diverging side wall in all three of the ducts they tested, at the lowest Reynolds number 
condition. No separation was observed at higher Reynolds numbers. Boundary layer 
measurements along the diverging side walls indicated the boundary layers were 
laminar at the lowest Reynolds number condition and turbulent at higher Reynolds 
numbers, explaining the existence of the separated flow region at the low Reynolds 
number. 

Flow with inlet swirl 

The efl^ect of inlet swirl on the transidon duct flow field is complex. On one hand, 
there appears to be little influence of inlet swirl on the duct surface static pressure 
distribution, shown in Figure V.39. However, the difference in detailed data, in 
particular the transverse components of velocity in all measurement planes for flow 
without and with inlet swirl is striking. For example, there is no evidence in plane 
3 A, Figure V.44, or plane 4A, Figure V.45 of the two pair of counter-rotating side 

wall vortices that were observed for flow without inlet swirl. 

* 
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When viewed looking upstream, the incoming flow with swirl is rotating clock- 
wise, resulting in a region outside the wall boundary layer of nearly constant negative 
axial vorticity. The overwhelming fluid flow effect was that the duct geometry driven 
convergence in the xz-plane and the divergence in the xy-plane was aided by the 
clockwise swirl in the upper right and lower left quadrants and opposed by the clock- 
wise swirl in the upper left and lower right quadrants, as shown in the sketch in 
Figure VI.9. This flow pattern established the static pressure distribution for plane 
1A seen in Figure V.46. The static pressure distribution again was saddle shaped, 
but the regions of highest static pressure were now located near the lower right and 
upper left comers and the regions of lowest pressure were located near the lower left 
and upper right comers. The static pressure gradient associated with the maximum 
and minimum static pressures drove the boundary layer flow near the lower right and 
upper left comers in two directions, some against the direction of rotation and some 



Figure VI.9 Sketch of swirl driven and duct geometry 
driven velocities for transition duct flow with inlet swirl 
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towards the direction of rotation. Near the lower left and upper right comers there 
was no radial static pressure gradient to balance the centrifugal forces and there the 
flow field was directed outward. Both of these effects can be seen in the plot of 
transverse velocity at plane 1A, Figure V.42. 

Hie response to the pressure gradient was again greater in the boundary layer 
where the momentum is less, resulting in skew-induced secondary flows. Near the 
four comers this produces positive vorticity. This can be seen in Figure VI. 10, which 
shows axial vorticity and transverse velocity components at plane 1A for flow with 
inlet swirl. The positive axial vorticity in the upper left quadrant resulted from the 
pressure gradi^ tangent ^ioihessuiiaee. Thepositiveaxial vorticity intheupper left 
quadrant resulted from the lack of radial pressure gradient. 

The saddle-shaped pressure distribution observed in plane 1A persisted to plane 
2A, shown in Figure V.47, however, it was not as pronounced. The complete 90° 
change in the orientation of the static pressure distribution seen in the flow field 
without inlet swirl was not observed. Rather, the regions of highest pressure remained 
in the lower right and upper left quadrants, but they appeared to be displaced counter 
clockwise, the highest pressure being nearer the y-axes and the lowest static pressures 
nearer the 2 -axes. 

Figure VI. 11 shows the axial vorticity at plane 2A. The positive axial vorticity 
observed in the upper left quadrant of plane 1A convected counter clockwise, against 
the principal direction of rotation, outwards towards the diverging side wall and out of 
the measurement region. This can be discerned from the surface oil film visualization 
shown In Figure V.40. The positive axial vorticity seen in the upper right quadrant 
of plane 1A has also convected outwards towards the diverging side wall, but this 
motion is clockwise, with the principal direction of rotation. 

At plane 3A the static pressure distribution, Figure V.48, does not differ sig- 
nificantly from the nonswirling static pressure at the same plane. The location of 
the regions of maximum static pressure were on the y-axes and the regions of mini- 
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mum static pressure were located on the z-axes. The static pressure distribution has 
changed its orientation nearly 90° from plane 1 A. The reversal of the static pres- 
sure distribution orientation for flow without inlet swirl was described earlier. For 
the nonswirling case, the reversal was accomplished by passing through a nearly flat 
pressure distribution. For the swirling case, the reversal of the pressure distribution 
orientation was accomplished without passing through a flat pressure distribution, 
the orientation of the static pressure distribution appeared to experience in effect a 
counter clockwise rotation. 

The axial vorticity at plane 3A is shown in Figure VI. 12. The positive axial 
vorticity in upper right quadrant was aided by the gradient in the static pressure 
distribution. The positive axial vorticity in the upper left quadrant was still out of 
the measurement region. This axial vorticity was not affected by the gradient in the 
pressure distribution because of its location. 

The static pressure distribution measured at plane 4A, Figure V.49, has the same 
shape as the distribution in plane 3A, however, the gradients were not as great. At 
this location and further downstream, the static pressure distribution had little effect 
on the boundary layer flow. Figure VI. 13 shows positive axial vorticity present in 
each comer at plane 4A. This was the same positive axial vorticity generated in 
each respective quadrant at plane 1A. The vorticity in the upper right quadrant was 
observed at all four measurement planes. The vorticity in the upper left quadrant was 
seen at plane 1A, but had convected out of the measurement region at planes 2A and 
3A, and reappears at plane 4A. 

The distortion of total pressure is most apparent in the lower right quadrant in 
planes 3 A and 4A, Figures V.48 and V.49. As in the nonswirling case, this distortion 
resulted from the convection of boundary layer fluid by secondary flows. The slight 
depression in total pressure near the centerline that appeared in all four planes is an 
artifact of the swirl generator which was described in Chapter V. All measurements 
indicate that the flow remains attached throughout the transition duct. 
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CHAPTER VII. CONCLUSIONS 

The combination of stationary blades followed by a rotating pipe containing a 
honeycomb and screen proved to be an effective way of generating a solid body 
swirling flow. The flow produced by the swirl generator was suitable for detailed 
downstream flow field measurements. The procedure of fabricating the stationary 
blades from cylindrical tubing with diminishing chord length saved considerable cost 
when compared to alternative fabrication methods involving extensive machine work. 

When utilized as exhaust system components of aircraft with rectangular nozzles, 
circular-to-rectangular transition ducts often involve an incoming flow that includes 
a swirling velocity component remaining from the gas turbine engine. Using the 
swirl generator improves the simulation of such transition duct flows. Detailed 
measurements of the transition duct flow field were obtained for flows with and 
without inlet swirl. A refinement of the five-hole probe aerodynamic measurement 
calibration and data reduction technique was used extensively for transition duct flow 
field measurements. These results provided greater insight into the fluid dynamics of 
circular-to-rectangular transition duct flows. 

Experimental results have shown that inlet swirl significantly changes the tran- 
sition duct flow field. Outside the boundary layer, the response of the flow field 
velocity to the changing duct geometry gives rise to the static pressure distribution. 
For nonswirling incoming flow the static pressure distribution produced skew-induced 
secondary flows within the boundary layer which evolved into two pair of counter- 
rotating side wall vortices. For the inlet swirl case the side wall vortices were not 
observed. The static pressure distribution was altered by the swirling velocity flow 
field to an extent that the side wall vortices were suppressed. The effects of inlet 
swirl should be included in the design of circular-to-rectangular transition ducts for 
aircraft exhaust systems. 

The ethylene trace gas technique provided an accurate Lagrangian flow field rep- 
resentation at high flow velocities within the transition duct. The methods developed 
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to obtain the ethylene trace gas data can be used in other experimental configurations 
with very little modification. In fact, this technique is presently being used in super- 
sonic flows to study mixing nozzles. The statistical methods used to help interpret the 
trace gas data provided a reliable way to quantify and compare contour characteristics 
such as shape, orientation or the amount of contour spreading. 
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CHAPTER VIII. RECOMMENDATIONS 
FOR FUTURE RESEARCH 

A requirement exists for further measurements in circular-to-rectangular transition 
ducts. Heat transfer is a particular problem in transition ducts used as components of 
aircraft exhaust systems. The ability of current CFD methods to predict heat transfer 
in transition ducts has to some degree suffered because of a lack of experimental 
heat transfer type data in these configurations. Experiments that utilize surface hot 
film gauges or liquid crystal methods to acquire surface shear stresses or convective 
heat transfer coefficients would help fulfill this requirement Inlet swirl should be 
included in this type of study. 

Another concern in aircraft propulsion design is the integration of various subsys- 
tems. A requirement exists to accurately predict how the performance of a propulsion 
component is influenced by the incoming flow resulting from upstream components. 
Using a swirl generator in the present transition duct study was partly a response to 
this concern. A continuation of transition duct research involves making flow field 
and surface measurements in a converging/diverging nozzle immediately downstream 
of a circular-to-rectangular transition duct. 

Another duct geometry that is interesting because of its application in aircraft 
propulsion systems is one with an S-shaped centerline. S-ducts are employed in both 
aircraft inlet and exhaust systems. A requirement exits for an S-duct experiment 
comparable to the present transition duct study. Swirl would again be a relevant inlet 
condition for an S-duct experiment 

Additional experiments and analysis should be directed towards the trace gas 
technique. There remains several unanswered questions regarding the interpretation 
of trace gas data when used to infer other flow field information. Using the trace 
gas technique to predict the effectiveness of transpiration and film cooling efforts in 
transition ducts and to provide a comparison with CFD predictions of these cooling 
methods should be explored. 
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APPENDIX A. A PROCEDURE FOR FIVE-HOLE 
PROBE CALIBRATION AND DATA REDUCTION 

The Five-Hole Probe 

Figure IV.2, reproduced below in Figure A.1, shows a typical five-hole probe and 
the nomenclature that will be used in this appendix. The openings of the four outer 
tubes were inclined, usually 45°, to their axis. The opening of the center tube was 
normal to its axis. This arrangement produced symmetry at the probe openings on 
the pitch and yaw planes. Away from the probe openings the five tubes were bundled 
inside the probe shaft. The tube openings were aligned with the probe shaft axis 
so the openings would not translate when the probe shaft is rotated. The pressure 
measured in the ith tube was denoted p,-. 




Figure A.l A typical five-hole probe 
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There are two methods used to obtain aerodynamic data using the five-hole probe, 
the yaw-nulling method and the stationary probe method. In the yaw-nulling method 
the probe was rotated about the probe shaft until the pressures p 2 and p 4 are equal. 
The yaw angle was measured directly. This appendix describes the calibration and 
data reduction procedure for the yaw-nulling method. 

The objective of five-hole probe calibration was to empirically determine the 
function, Equation A.l, relating the flow total pressure po* the flow static pressure p, 
and pitch angle a, to the subspace (p,-,pj,pjt) of the five measured pressures. 

F : (p 0 ,p, a) -+ (pi,Pj,Pk) C (pi, . . . ,P5) (A.l) 

The inverse of this function, A.2 , used in data reduction, allowed the total and static 
pressure and pitch angle to be calculated from three of the five measured pressures. 

: (Pi,Pj,Pk) -* (P0.P,a) (A.2) 

The success of this procedure depended on the existence of this inverse function. 
The conditions necessary for the existence of the inverse function are described in 
the following section. 

Calibration 

Experimental procedure 

The calibration was accomplished by placing the probe in a uniform flow and 
varying the probe pitch angle and Mach number while measuring the probe’s five 
pressures. By varying the Mach number the ratio of static to total pressure was 
varied. This procedure was performed in the Internal Fluid Mechanics Facility. The 
probe and actuator were mounted on a gimbal with the probe placed so the probe 
openings were located at the gimbal ’s center of rotation. This allowed the pitch 
angle to change without a translation of the probe openings. While the kth Mach 
number remained constant the probe was positioned at the Zth different pitch angle 
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and the five probe pressures, probe pitch angle, and the flow total pressure and static 
pressure were recorded. Measurements were recorded at a total of AT* different pitch 
angles for the kth Mach number. This procedure was repeated for a total of ac h 
different Mach numbers. 

Dimensional considerations 

The five pressures measured by the probe were nondimensionalized by Equation 
A.3. This is the same pressure coefficient given by Equations V.l and V.2 that were 
used to nondimensionalize experimental static and total pressure data. 

= (A.3) 

PO-P \ d) 

The pressure coefficients could depend on flow parameters such as the Mach number 
M, the Reynolds number Re, the specific heat ratio 7, and probe parameters such 
as the pitch angle a, and the probe length-to-diameter ratio l/d. Air was used for 
both the calibration and experiments so the specific heat ratio was removed from the 
list of independent variables. The Mach and Reynolds number could not be varied 
independently in the Internal Fluid Mechanics Facility so Mach number dependence 
was considered while the Reynolds number was removed from the list of independent 
variables. 

When calibrating the five-hole piobe it was the pitch angle relative to the probe 
shaft a, that was varied and recorded, but it was the pitch angle relative to the probe 
openings 0, to which the probe responded. These different angles are shown in Figure 
A.2. The difference between the two pitch angles was the total pitch difference angle 
0o, where 80 = 8 - a. Figure A.3 shows typical calibration data results. The five 
pressure coefficients are plotted as functions of 8 for three different Mach numbers. 
These results are essentially independent of Mach number for the range of Mach 
number tested. It was in the angle 0o that the Mach number and length-to-diameter 
ratio dependence appeared, as symbolized by Equation A.4. 


Si = f>( 8 ) = /,( a + 00),* = 1 , ... ,5 00 = 00 (m, 2 


(A.4) 



160 



Two principal causes contributed to the total pitch difference angle Qq. The first 
was any difference from perpendicular in the angle between the probe sensing area 
and the probe shaft caused by probe construction. This difference in angle could be an 
intentional design feature of the probe or could result unintentionally from imperfect 
construction. This was the static pitch difference angle. The static pitch difference 
angle was unaffected by flow Mach number and the probe length-to-diameter ratio. 
The second cause was the deflection of the probe produced by aerodynamic loading, 
which did depend on flow Mach number and the probe length-to-diameter ratio. This 
was the dynamic pitch difference angle. The sum of these two angles was the total 
pitch difference angle 0q. 

The pitch angle 0, not a, was used for probe calibration and data reduction. Using 
the pitch angle 6 had several important advantages. 

1. The calibration results were independent of Mach number for a significant range 
of subsonic Mach numbers, as seen in Figure A.3. 
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pitch ongle 8 


Figure A.3 Nondimensional pressure functions for a typical five-hole probe 

2. The symmetry of the probe lead to several important simplifications in the analysis 
of the calibration data and the data reduction procedure. This is discussed in 
detail later. 

3. During experimental tests probes are occasionally bent slightly by striking solid 
surfaces. When this occurs, experience has shown the functions A.3 are unaf- 
fected, only the total pitch difference angle, Bo, is changed. Recalibrating the 
probe was often unnecessary. 

The objective of probe calibration, to relate po, p, and a to three of the five 
measured pressures remained unchanged, however, the functional relation actually 
used in the calibration and data reduction procedure is more accurately represented 
by Equations A.5 and A.6. In order to satisfy the original objective the angle 6q must 
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be known to determine a. The procedures used to find 8q for calibration and data 
reduction are described later. 

F : (PQ,P,0) -*• ( Pi,Pj,Pk ) C (pi, . . . ,P5) (A.5) 

F~ l : (pi,Pj,Pk) -+ (pO,P,0) (A.6) 

Symmetry conditions and Taylor’s series 

The symmetry of the probe required the pressure coefficients satisfy the conditions 
given in Equations A.7 to A.ll. 

/l(«) = M-f) (A.7) 


h«>) = h(-0) (A.8) 

M«) = M«) (A.9) 

M») = /,(-«) (A.10) 

fi(0) = MO) = MO) = MO) (A.ll) 

For subsonic flow, aerodynamic considerations required Equation A.12 be satisfied. 

MO) = 1 (A.12) 


The Taylor’s series for the pressure coefficients is shown in Equation A. 13. 

OO 

fi — a i,Q + Oj',1^ + O-i, 2^ 2 + • ‘ ‘ = a i,j& ■> i = 1, . . . , 5 

i= i 


(A. 13) 
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The symmetry conditions had important consequences for the Taylor’s series of each 
pressure coefficient. Condition A.7 required Equations A.14 and A.15 be satisfied. 

o,\ t j — Ajj) j — 0 ? 2, . . . (A.14) 

o\ ,j = J = 1,3,... (A.15) 

Condition A. 8 required Equation A. 16 be satisfied. 

0 - 2 , j = 0, j = 1, 3, . . . (A.16) 

Condition A.9 required Equation A. 17 be satisfied. 

a 2 , j = 04 , j, j = 0,1,... (A. 17) 

Condition A. 10 required Equation A. 18 be satisfied. 

Oh,j = 0, j = 1, 3, . . . (A.18) 

Condition A. 11 required Equation A. 19 be satisfied. 

ai,i = a«, 2 = = ai, 4 (A. 19) 

The aerodynamic condition, A. 12, required Equation A.20 be satisfied. 

05,0 = 1 (A.20) 

Using these results, the Taylor’s series of the pressure coefficients are given in 
Equations A.21 through A.25. 

fi(0) = ai,o + ai,!# + ai,20 2 H (A.21) 


/ 2 ( 0 ) = 01,0 + a 2 ,20 2 + • • • 


(A. 22 ) 
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h(Q) = ai,o - + <*i,20 2 H — 


(A.23) 


/4 (0) = <*1,0 + <*2,20 2 + • • • 

(A.24) 

/s(0) = 1 + a5,20 2 + • • • 

(A.25) 


Inverse function theorem 


For the inverse function given by Equation A.6 to exist the inverse function 
theorem required three of the five measured pressure satisfy Equation A.26 for all 
values of po, p, and 6 of interest. 



(A.26) 


The magnitude of this determinant was an indication of the suitability of the 
inverse function for data reduction. A large magnitude was desirable. Substituting 
the pressure coefficients into Equation A.26 resulted in Equation A.27. The case when 
po = p corresponded to no fluid motion and was of no importance. By using the 
Taylor’s series of the pressure coefficients in Equation A.27 it was possible to predict 
which three pressures are most suitable for use in calibration and data reduction. 


dp, dp, 


det 


d££ 

dpo dp d6 

& 

apt) dp ad 

dp* dpi t dp* 

ciphi Up ~W 


= (PO “ P) 


dfi 

de 


dfk 


dfj. 


L U ~h)~ ft - fi) + -£{fk - /.) 


de 


(A.27) 


The magnitude of the determinant in Equation A.27 for the ten combinations 
resulting from five pressures taken three at a time is summarized in Table A.l. Three 
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Table A.l The ten pressure combinations and the 
corresponding value of the determinant given in Equation A.27 


Pressure Combinations 

' dp t dpi dpt ' 

dpo dp of 

** & % & 

dp* dp* dpk 

- dpo 5p a<T - 

P 1 ,P 2 ,P 3 

(po - p)|2ai ) i(a2,2 - ai,2)<? 2 | 

P 1 ,P 2 ,P 4 

0 

P 1 ,P 2 ,P 5 

(PO -p)| ~ ai,oai,l+ 
201,0(02,2 ~ «l, 2 ) 0 + 
a 1,1(02,2-05,2)^1 

P 1 ,P 3 ,P 4 

(po-p)|2ai,i(a 2 ,2 - oi,2)^ 2 | 

Pl,P 3 ,P 5 

(PO -p)| - 2 ai,oai,i+ 
2 ai,i(ai ,2 - o 5 , 2 ) 0 2 | 

P 1 ,P 4 ,P 5 

(PO - p) |-ai,oai,i + 

201,0(02,2 - 01,2)^+ 
01,1(02,2 - O,«j,2)0 2 | 

P 2 ,P 3 ,P 4 

0 

P 2 ,P 3 ,P 5 

(PO -p)| - 01,001,1+ 
2oi,o(a2,2 - 01,2)6!+ 

ai, 1(02,2 - O5,2)0 2 | 

P 2 ,P 4 ,P 5 

0 

P 3 ,P 4 ,P 5 

(po -p)| - ai, 0 oi,i + 
201,0(02,2 - 01,2)6!+ 
01,1(02,2 - a5,2)^ 2 | 


combinations, pi,P2,P4» P 2 ,P 3 ,P 4 , and p2,P4,P5, produced a determinant that was 
zero for all values of 6 . The determinant of two combinations, pi, P 2 , pa and pi , pa, P4 
vanished when 0 = 0. The remaining five combinations all had nonzero determinants 
in the neighborhood of 8 = 0. Of these, the combination pi,pa,Pf, was the most 
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suitable to use for calibration. The magnitude of the determinant for this combination 
was twice as large as the remaining four combinations at 6 = 0. 

Determining 9o 

To calculate the angle 9 for each calibration measurement it was necessary to 
determine the total pitch difference angle 0o,* for each Jfcth Mach number at which 
calibration measurements were recorded. This was done numerically by calculating 
the value of a which satisfies f^ k == /a,*. This was accomplished by finding the 
least squares linear approximations given by Equation A.28. 

h,k,i ~ h,k,i ~ co,* + ci,*a (A.28) 

The coefficients co,* and ci,* were calculated by solving the system of Equations 
A.29 and A.30. In these equations /i,*,/, /a,*,/ and a k j are determined from the 7 th 
measurements recorded at the kth Mach number. 

N k N k 

CQ,kN k + Cl,* Y <*k,i = Y Cfr.M “ h,k,i) (A.29) 

/=i /=i 

N„ N k 

CO,* Y a k,l + C 1 ,k Y a k,l = (fhk,l - h,k,l)<*k,l (A.30) 

1=1 1=1 

Oo k was given by Equation A.31. 

Cn Jb 

Oo, k = (A.31) 

ci,* 

Least squares approximation 

To approximate the three pressure coefficients /i, /a and /s with the first three 
terms of their Taylor’s series, only four of the nine coefficients must be determined; 
Qi,o» cr l , * , ai ,2 and 05 2 - This was done using a least squares procedure. The error 
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of the approximation, for all TV* experimental measurements recorded at all Aj^ach 
Mach numbers, was defined by Equation A.32. 

^Mach N k 2 

E = ^2 E [/i.M " ( ai >° + ai > 1 $k ,i + a htfk,i)] 

Jk=i t=i 
*Mach N k 

+ E E [h> k >l - ( ai -° “ a hrfk,l + ai,2 Ql,i)] (A.32) 

Jk=i 1=1 

^Mach Nk 

+ E E [hk,l~ (i+a^.i)] 

*=i 1=1 


The least squares procedure required finding the four coefficients which minimize 
this error. These four coefficients satisfy Equation A.33. The four coefficients were 
found by solving Equations A.34 through A.37. Equations A.35 and A.37 were solved 
independently, Equations A.34 and A.36 were solved simultaneously. 


dE _ dE_ _ dE_ _ dE 
da 1,0 dau dai,2 dag, 2 


(A.33) 


^Mach ^Mach N k ^Mach N k 

a,,» £ N k +a u £ £f>?,,= r £ £(/u,l + / W ) (A.34) 

k=l k= 1 1=1 Jfc=l 1=1 

A Mach w* 1 ^Mach N k 

a M E E^U = 5 E Eihv-hk^i (A.35) 

fc=l 1=1 fc=l 1=1 


ai,0 


^Mach N k 

E EX/+ a 1,2 


*=1 1=1 


^Mach iv* 

E E* 


jt=i 1=1 



^Mach jy* 

E E (fw + h,k,i) e h ( A -36) 


fc=i 1=1 


*Mach ^Mach iy* 

E E 6 k,i + as, 2 E E e k,i = 

*=1 1=1 k—l 1=1 


*Mach N k 

E ^2h,k,i e k,i 

k=l 1=1 


(A.37) 
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Data Reduction 

The data reduction procedure used the truncated Taylor’s series approximations 
of the pressure coefficients determined in the calibration procedure to calculate the 
flow field total pressure, static pressure, and direction from the probe pressures and 
yaw angle measured during an experiment. This required finding the inverse function, 
A.6, whose existence was established in the previous section. 

Equation A.38, a ratio of differences of measured pressures, was used to determine 
9. Any ratio of differences of measured pressures is equivalent to the correspond- 
ing ratio of differences of pressure coefficients, making the ratio of differences of 
measured pressures independent of total and static pressure and dependent on 9 only. 
For a ratio of differences of pressure coefficients, the denominator of each pressure 
coefficient, po - p, appears in both the numerator and denominator of the ratio of 
differences and cancel. The static pressure term which appears in the numerator of 
each pressure coefficient, p, - p, is cancelled by the difference operations in the ratio 
of differences. These cancellations leave only the measured pressures, pi, in the ratio 
of differences of pressure coefficients. 

^(pi — pa) = lifi ~ h) = <*i,i* (A<38) 

P5 - i(pi + pa) h - \(h + h) («6,2 - <*i,2 )0 2 + 1 - <*i,o 

Equation A.38 is a quadratic in 8, which appears contrary to the conclusion based 
on the inverse function theorem that the function A.5 was single valued. However, the 
function A.5 was single valued within the range of pitch angles used for calibration. 
Experience has shown that one root of Equation A.38 lies within the range of pitch 
angles while the second root lies well outside the range of pitch angles. In practice 
there has been no difficulty deciding which root of Equation A.38 is correct. 

The value of 6 found using A.38 was used in Equation A.39 to calculate the total 
and static pressure. 


<*1,0 + <*1,20 2 

1 - ai,o - ai ( 20 2 

po' 


1(P1+P3)' 

1 + a.5,20 2 

“<*5,2 0 2 

.P. 


P5 


(A.39) 
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These values were then used in Equation A.40 to determine the Mach number. 



(A.40) 


Figure A.4 shows a nulled five-hole probe and the flow direction relative to the 
probe coordinate system. In this coordinate system the z'-axis is parallel with the 
probe shaft. The z'-axis is parallel to the probe pitch plane when ip = ip$. The angle 
ipo is the yaw reference angle. The components of a unit length vector in the flow 
direction relative to the probe coordinate system are given by Equation A.41. 

U X 1 

Uy' 

_U Z I 

To determine the flow direction relative to the probe coordinate system required 
knowing the total pitch difference angle do and the yaw reference angle ipo. Two 
ways of finding the angle Qo are: 

1. Make a reference flow measurement at a location where the flow direction is 
known relative to the probe coordinate system. 

2. Use the values of $o found during calibration to estimate the value of 9q in the 
survey. 

If a reference measurement was used to obtain $o the measurement should be 
made under conditions where the Mach number and probe length-to-diameter ratio 
are nearly equal to the conditions expected in the survey. Usually, this reference 
measurement is made at a location in the experiment where the flow direction is 
known a priori, or in a different facility where the flow direction is known. The pitch 
angle 9, calculated from the reference pressure measurements, and the known pitch 
angle a are used to calculate 0q. 

Estimating the angle 0q from calibration data was difficult because the angle 
$o can change if the probe has been bent slightly after calibration. The reference 


cos (ip — ipo) cos (6 — 0o) 
sin (ip — ipo) cos (0 — <?o) 
sin(0 — 0o ) 


(A.41) 
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Figure A.4 Flow direction relative to the probe coordinate system 
measurement method was unaffected by any changes in the angle 8q occurring prior 
to the reference measurement. 


The yaw reference angle -00 was necessary for practical considerations. Typically, 
yaw angles are measured from an arbitrary origin. When the reference measurement 
was recorded, 0o was assigned the value of the of yaw angle at which the probe 
nulled. This allowed all future yaw angles to be referred to this angle. 

To determine the flow direction relative to the experiment coordinate system 
required knowing the orientation of the probe relative to the experiment coordinate 
system. This orientation was determined by three angles, 4>\, <j> 2 , and fa, shown in 
Figure A.5. The first angle, fa, is a rotation about the 2 -axis, the second angle, fa, is 
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a rotation about the new y-axis, the third angle, <f> 3, is a rotation about the new x-axis. 
Equation A .42 was used to transform the flow direction from probe coordinates to 
experiment coordinates. 


u x 


Uy 

m U z m 

— 


cos^i 
sin <f>\ 

0 


1 

0 

0 


— sin^i 0 
cos 0 

0 1 

0 0 

cos <^3 — sin <f> 3 

sin ^3 cos ^3 


cos <f >2 0 

0 1 
sin <j> 2 0 



1 

1 


Uy, 




— sin <l> 2 
0 

cos <j>2 J 


(A. 42 ) 



Figure A .5 The relation between the probe coordinate 
system and the experiment coordinate system 
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A Fortran Program For Five-Hole Probe Calibration 


c I I I I I I l-C 

c c 

c PCAL 1.1 (03/09/90) c 

c c 

c Bruce Reichert c 

c Mail Stop 5-11 c 

c NASA Lewis Research Center c 

c 21000 Brookpark Road c 

c Cleveland OH 44135 c 

c (216) 433-8397 c 

c c 

c I [ | purpose | I I — 1 -c 

c c 

c The purpose of this program is to determine the four coefficients c 

c required to approximate the nondimensional pressures measured c 

c during the five-hole probe calibration procedure . c 

c c 

c I I Input - Output I i | -c 

c c 

c The input file is the calibration data. This is an ASCII file c 


c named WW_XX .MM_DD_YY where the WW XX is the probe serial number c 
c and MM_DD_YY is the date of the calibration. The first record is c 
c an integer equal to the number of different Mach numbers for which c 
c calibration data were recorded. The calibration data follows c 
c grouped by Mach number. The first record is the number of different c 
c pitch angles recorded at the particular Mach number followed by the c 


c Mach number. This is followed by the data recorded for each pitch c 

c angle. The format for this data is (2 (2x, f6.2), 7 (2x, f7.3)) and c 

c is written in the order alpha, psi, pi, p2, p3, p4, p5, p0, p. c 

c c 

c The ouput file contains the four calibration coefficients. Their c 

c format is (4f) and are written in the order alO, all, al2, a52. c 

c c 

c — I j Variables i I I -c 

c c 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c- 


Real c 

c 

alO, all, al2, a52 - The coefficients of the third order c 

approximation of the nondimensional pressures as c 

functions of theta. c 

alpha - The pitch angle relative to the probe shaft. c 

amat - A matrix used to calculate least square coefficients. c 

bvec - A vector used to calculate least square coefficients. c 

cO, cl - The coefficients used to calculate thetaO . c 

mach - An array containing each flow field Mach number. c 

p - Flow field static pressure. c 

pO - Flow field total pressure. c 

pi, p2, p3, p4, p5 - The five measured probe pressures. c 

theta - The pitch angle relative to the probe openings. c 

thetaO - An array containing the total pitch difference angle for c 
each Mach number. c 

c 

Integer. c 


c 

n - An array containing the number of different pitch angles c 


recorded for each Mach number. c 

nmach - The number of different Mach numbers. c 

c 

Character c 

c 

date - The character representation of the five-hole probe c 

calibration date. c 

infile - The name of the input file c 

outfil - The name of the output file c 

serial - The character representation of the five-hole probe serial c 

number. c 

c 


f -c 


real amat(4,4), bvec(4), thetaO (10), mach(10) 


integer n(10) 
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character*5 serial*5, date*8 
character*80 infile, outfil 

gamma = 1.40 

write (*, ' (a)') ' Enter probe serial number' 
read(*, '(a)') serial 

write (*, ' (a ) 9 ) ' Enter probe calibration date' 
read ( *, ' (a) ' ) date 

infile = ' /homel/f sreich/calibration/original/' 

1 //serial//' .' //date 

outfil = ' /homel/f sreich/calibration/parameter/ ' 

1 //serial//' .' //date 

open (unit = 1, file = infile, status - 'old') 
read (1, ' (i) ' ) nmach 

c For each Mach number find the difference between the 
c' pitch angle relative to the probe shaft, alpha, and the pitch angle 
c relative to the probe openings, theta. This is the total pitch 
c difference angle, thetaO. 

do 10 k = 1, nmach 

read {1, ' (i) ' ) n (k) 
read(l, '(f)') mach(k) 

do 20 1 = 1, n (k) 

read(l, '(2(2x, f6.2>, 7 (2x, f7.3))') alpha, psi, pi, p2, P 3, 
1 p4, p5, p0, p 

fl = (pi - p) / (p0 - p) 
f 3 = (p3 - p) / (pO - p) 
amatd, 1) = amat(l, 1) + alpha 

amatd, 2) = amatd, 2) + 1.0 

amatd, 1) = amat(2, 1) + alpha ** 2.0 
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amat(2, 2) = amat(l, 1) 
bvec(l) = bvec(l) + (fl - f3) 

bvec (2) = bvec (2) + (fl - f3) * alpha 

20 continue 

det = amat(l, 1) * amat(2, 2) - amat(2 r 1) * amat(l, 2) 

cO = (amat (1,1) * bvec(2) - amat (2, 1) * bvec(l)) / det 

cl = (bvec ( 1 ) * amat ( 2 , 2 ) - bvec (2) * amat(l, 2)) / det 

thetaO(k) = cO / cl 

amat (1,1)=0.0 
amat (1, 2) = 0.0 
amat (2, 1 ) =0.0 

amat (2, 2) =0.0 
bvec ( 1 ) =0.0 

bvec (2) =0.0 

10 continue 

rewind (unit = 1) 

c Curve fit the three dimensionless pressures as a funtion 
c of true pitch angle. 

read(l, Mi)') nmach 
do 30 k = 1, nmach 

read (1, Mi) ' ) n (k) 
read(l, Mf)') mach(i) 

do 40 1 = 1, n (k) 

read(l, '<2(2x, f6.2), 7 (2x, f7.3))') alpha, psi, pi, P 2, P 3, 
1 p4, p5, pO, p 

fl = (pi - p) / ( P 0 - p) 

f 3 = (p3 - p) / (p0 - p) 

f 5 = (p5 - p) / (pO - p) 

theta = alpha + thetaO(k) 

amat ( 1 , 1) = amat(l, 1) + 1.0 
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amat (1, 

2) 

= amat (1 , 

amat (2, 

2) 

.= amat (2, 

bvec (1) 


- bvec(l) 

bvec (2) 


= bvec (2) 

bvec (3) 


= bvec (3) 

bvec (4 ) 


= bvec (4) 


40 continue 

30 continue 


2) + theta ** 2.0 
2) + theta ** 4.0 
+ 0.5* (fl + f3) 

+ 0.5 * (fl + f3) * theta ** 2.0 
+ 0.5 * (fl - f3) * theta 
+ f 5 * theta ** 2.0 


amat (2, 1) 
amat (3, 3) 
amat (4 , 4 ) 


= amat (1, 2) 
= amat (1, 2) 
= amat«(2, 2) 


det = amat(l, .1) * amat (2, 
alO = (bvec(l) * amat (2,2) 
all = bvec(3) / amat (3, 3) 
al2 = (amat (1,1) * bvec (2) 
a52 = (bvec(4) - amat (1,2)) 


2) - amat (2 , 1) * amat(l, 2) 

- bvec (2) * amatd, 2)) / det 

- amat {2 r 1) * bvec(l)) / det 
/ amat (4, 4) 


open (unit = 2, file = outfil, status = 'new') 
write {2, 9 ( 4 f ) ' ) alO, all, al2, a52 


stop 

end 


A Fortran Program For Five-Hole Probe Data Reduction 


c- 

c 

c 

c 

c 

c 


PDATA 1.1 (03/09/90) 

Bruce Reichert 
Mail Stop 5“ 1 1 


I -c 
c 
c 
c 
c 
c 
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c NASA Lewis Research Center c 

c 21000 Brookpark Road c 

c Cleveland OH 44135 c 

c (216) 433-8397 c 

c c 

c I j | purpose |— | ! I -c 

c c 

c The purpose of this program is to read the pressures and yaw angles c 
c measured by the five-hole probe and use the calibration data for c 
c that probe to calculate the total and static pressure, Mach number, c 
c and direction cosines for each fivehole probe reading. c 

c c 

c I I Input - Output j I I — c 

c c 

c The first input file is an ASCII file named WW_XX.MM_DD_YY. WW__XX c 
c is the five-hole probe serial number and MM_DD_YY is the date of c 
c probe calibration. This file contains the four coefficients, alO, c 
c all, al2, a52, written in the format (4f) . The second input file c 
c is an ASCII file that contains the experimental data. The first c 
c record, an integer, is the number of different locations at which c 
c data were recorded. This is followed by the data for each location. c 
c The format for this data is (2x, f 7 .3, 2x, f 5 . 1, 5 (2x, f7 . 3) ) and is c 
c written in the order h, psi, pi, p2, p3, p4, p5 . The program c 

c expects the first data location to be the reference reading used c 
c to find the reference pitch angle, thetaO, and the reference yaw c 
c angle, psiO. c 

c c 

c The output file is an ASCII file that contains the reduced data. c 
c The first record, an integer, is the number of locations at which c 
c data were reduced. This is one less than the number of data c 

c locations in the second input file. This is because the reference c 
c reading is not a data location for data reduction, it only used as c 
c a reference for other data locations. The integer is followed by c 
c the data for each location written in the order x, y, z, pO, p, c 

c mach, ux, uy, uz . The format is (9 (2x, f 7 . 3) ) . c 

c c 

c | j Variables I I | -c 

c c 
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c Real c 
c c 
c amat - The matrix used to calculate total and static pressure. c 
c bvec - The vector used to calculate total and static pressure. c 
c gamma - The specific heat ratio. c 
c h - The distance from the probe to the origin of the probe c 
c coordinate system. c 
c mach - The Mach number. c 
c p - The static pressure. c 
c phil, phi2, phi3 - The three angles which determine the orientation c 
c of the probe coordinate system with respect to the c 
c experiment coordinate system. c 
c pO -* The total pressure. c 
c ux,uy,uz - The direction cosines of the flow direction with respect c 
c to the experiment coordinate system. c 
c uxl,uyl,uzl - Intermediate results used in transforming the c 
c ux2,uy2,uz2 direction cosines from the probe coordinate system c 
c . to the experiment coordinate system. c 
c ux3,uy3,uz3 - The direction cosines of the flow direction with c 
c respect to the probe coordinate system. c 


c x,y,z - The cartesian coordinates with respect to the experimental c 


c coordinate system of the location at which the data were c 
c recorded. c 
c xO,yO,zO - The location of the origin of the probe coordinate c 
c system with respect to the experiment coordinate system. c 
c xl,yl,zl - Intermediate results used in transforming the. data c 
c x2,y2,z2 location coordinates from the probe coordinate system c 
c to the experiment coordinate system. c 
c x3,y3,z3 - The cartesian coordinates with respect to the probe c 
c coordinate system of the location at which the data were c 
c recorded . c 
c c 
c Integer c 
c c 
c ndata - The number of locations at which data were recorded. c 
c plane - The number of the data plane at which the data were c 
c recorded . c 
c slot - The number of the probe opening at which the data were c 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c- 


recorded. c 

c 

Character c 

c 

date - The character representation of the five-hole probe c 

calibration date. c 

machno - The character representation of the freestream Mach number c 
at which data were taken. c 

outfil - The name of the output file. c 

planno - The character representation of the number of the data c 

plane at which the data were recorded. c 

serial - The character representation of the five-hole probe serial c 
number. c 

slotno - The character representation of the number of the probe c 
opening at which the data were recorded. c 

c 

| j | j | I I — c 


real amat ( 2 , 2 ) , bvec ( 2 ) , mach 


integer plane, slot 


character serial*5, date*8, machno*7, planno*7, slotno*6 
character*80 infile (2), outfil 

pi = atan2 (0.0, -1.0) 
gamma - 1.40 

write (*, ' (a) ' ) ' Enter probe serial number' 
read ( * , ' (a) ' ) serial 


write(*, '(a)') ' Enter probe calibration date' 
read (*, ' (a) ' ) date 

write (*, '(a)') ' Enter Mach number' 
read (* , '(f)') mach 

if (mach .eq. 0.20) machno = 'mach20/' 

if (mach .eq. 0.35) machno = 'mach35/' 

if (mach .eq. 0.50) machno = 'mach50/' 
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write ( * , '(a)') ' Enter plane number' 
read (*, ' (i) ' ) plane 

if (plane .eq. 1) planno = 'planel/' 

if (plane .eq. 2) planno = 'plane2/' 

if (plane .eq. 3) planno = 'plane3/' 

if (plane . eq. 4) planno = 'plane4/' 


write(*, ' (a)M ' Enter slot number' 
read ( *, ' (i) ') slot 


if 

(slot 

.eq. 

1) 

slotno 


' slotOl' 

if 

(slot 

.eq. 

2) 

slotno 


' slot02' 

if 

(slot 

.eq. 

3) 

slotno 

= 

' slot03' 

if 

(slot 

.eq. 

4) 

slotno 

= 

' slot04 ' 

if 

(slot 

.eq. 

5) 

slotno 

= 

' slot05' 

if 

(slot 

.eq. 

6) 

slotno 

= 

' slot06' 

if 

(slot 

.eq. 

7) 

slotno 

= 

' slot07' 

if 

(slot 

.eq. 

8) 

slotno 

= 

' slot08 ' 

if 

(slot 

.eq. 

9) 

slotno 

= 

' slot09' 

if 

(slot 

.eq. 

1°) 

slotno 

= 

' slot 10 ' 

if 

(slot 

.eq. 

11) 

slotno 

= 

' slot 11' 

if 

(slot 

.eq. 

12) 

slotno 

= 

' slot 12 ' 

if 

(slot 

.eq. 

13) 

slotno 

= 

' slotl3' 

if 

(slot 

.eq. 

14) 

slotno 

= 

' slot 14 ' 

if 

(slot 

.eq. 

15) 

slotno 

= 

' slot 15 ' 


infile ( 1 ) 

1 

infile (2) 

1 

outf il 

1 


' /homel/f sreich/calibration/parameter/' // serial 
// ' .' // date 

' /homel/f sreich/pressur e/ original/' // machno 
// planno // slotno 

' /homel/fsreich/pressure/reduced/' // machno 
// planno // slotno 


open (unit = 1, file = infile(l), status = 'old') 
open (unit = 2, file = infile (2) , status = 'old') 
open (unit = 3, file = outfil, status = 'new') 


write ( * 
read ( * , 


'(a)') ' Enter x reference location, xO 
(f)') xO 
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write {*, '(a)') ' Enter y reference location, yO' 
read (* , ' (f ) ' ) yO 

write (*, '(a)') ' Enter z reference location, zO' 
read (*, ' (f ) ' ) zO 

write(*, '(a)') ' Enter. first transformation angle, phil' 
read ( * , '(f)') phil 

write (*, '(a)') ' Enter second transformation angle, phi2' 
read ( *, '(f)' ) phi2 

write(*, ' (a) f ) ' Enter third transformation angle, phi3' 
read {* , '(f)') phi3 

c Read the calibration coefficients. 

read (1, ' <4f) ' ) alO, all, al2, a52 

read (2, ' (i ) ' ) ndata 
write(3, '(i)') ndata - 1 

c The first reading is expected to be the reference reading. 

read (2, ' (2x, f7 . 3 , 2x, f 5 . 1 , 7 (2x, f 7 . 3) ) ' ) h, psiO, pi, p2, p3, p4, 
1 p5, ptotal, pstatic 

c Calculate the reference pitch angle, thetaO, and the reference 
c yaw angle, psiO. 

g = 0.5 * (pi - p3) / (p5 - 0.5 * (pi + p3 ) ) 
a * g * (a52 - al2) 
b = - all 

c = g * (1.0 - alO) 
if (abs (g) .It. 1 . 0e-008 ) then 
thetaO =0.0 
else 

thetal = {- b + sqrt (b * b - 4.0 * a * c) ) / (2.0 * a) 
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theta2 = (- b - sqrt <b *b-4.0*a*c)) / (2.0 * a) 
if (abs(thetal) .It. abs(theta2)) then 
thetaO = thetal 
else 

thetaO = theta2 
end if 
end if 

thetaO = thetaO * pi / 180.0 
psiO = -psiO * pi / 180.0 

do 20 i = 1, ndata - 1 

read (2, ' (2x, f7 .3, 2x, f 5. 1, 5 (2x, f7 . 3) ) ' ) h, psi , pi, p2, p3, p4 , 
c p5 

c Calculate the pitch angle, theta. 

g = 0.5 * (pi - p3) / (p5 - 0.5 * (pi 4 - p3) ) 
a = g * (a52 - al2) 
b = - all 

c=g*(1.0-a!0) 


if (abs (g) . It . 1 

. 0e-008 ) 

then 






theta = 0.0 








else 








thetal = (- b + 

sqrt (b 

* b - 4.0 

* a * c) ) 

/ 

(2. 

.0 * 

a) 

theta2 = (- b - 

sqrt (b 

* b - 4.0 

* a * c) ) 

/ 

(2. 

.0 * 

a) 


if (abs(thetal) .It. abs(theta2)) then 
theta = thetal 
else 

theta = theta2 
end if 
end if 

c Use the pitch angle, theta, to determine the coefficients of the 
c linear equation used to determine the total and statid pressure. 


amat(l, 1) = alO + a!2 * theta ** 2.0 
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A 


amat ( 1 , 
amat (2, 
amat (2 f 
bvec (1) 
bvec (2) 


2) = 1.0 - alO - al2 * theta ** 2.0 

1) = 1.0 + a52 * theta ** 2.0 

2) = -a52 * theta ** 2.0 
= 0.5 * (pi + p3 ) 

= p5 


c Calculate the total and static pressure by solving the linear 
c equation. 


det - amat(l, 1) * amat (2, 2) - amat (2, 1) * amat(l, 2) 
p0 = (bvec(l) * amat (2, 2) - bvec (2) * amat(l, 2)) / det 
p = (amat(l f 1) * bvec (2) - amat (2, 1) * bvec(l)) / det 

c Calculate the Mach number from the total and static pressure. 


mach = sqrt ( (2.0 / (gamma - 1.0)) * (((p0/p) ** ((gamma - 1.0) 
1 / gamma) ) - 1 . 0) ) 

c Convert the pitch and yaw angles from degrees to radians. 


theta = theta * pi / 180.0 
psi = -psi * pi / 180.0 

c Calculate the direction cosines and data location relative to 
c the probe coordinate system. 

ux3 = cos (psi - psiO) * cos (theta - thetaO) 
uy3 = sin (psi - psiO) * cos (theta - thetaO) 
uz3 = sin(theta - thetaO) 

x3 = 0 
y3 =0 
z3 = h 

c Transform the direction cosines and data location from the probe 
c coordinate system to the experiment coordinate system. This is 
c done in three transormations corresponding to the three angles, 
c phil, phi2 , phi 3, which define the orientation of the probe 
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c coordinate system relative to the experiment coordinate system. 
ux2 = ux3 

uy2 = uy3 * cos(phi3) - uz3 * sin(phi3) 
uz2 = uy3 * sin(phi3) + uz3 * cos(phi3) 

x2 = x3 

y2 = y3 * cos(phi3) - z3 * sin(phi3) 
z2 = y3 * sin(phi3) + z3 * cos(phi3) 

uxl = ux2 * cos (phi2 ) - uz2 * sin(phi2) 
uyl = uy2 

uzl = ux2 * sin(phi2) + uz2 * cos(phi2) 

xl = x2 * cos(phi2) - z2 * sin(phi2) 
yl - y2 

zl = x2 * sin(phi2) + z2 * cos(phi2) 

ux = uxl * cos(phil) - uyl * sin(phil) 
uy = uxl * sin(phil) + uyl * cos(phil) 
uz = uzl 

x = xO + xl * cos(phil) - yl * sin(phil) 
y = yO + xl * sin(phil) + yl * cos(phil) 
z = zO + zl 

write (3 , ' (9 (2x, f7 .3) ) ' ) x, y, z, pO, p, mach, ux, uy, uz 
20 continue 


stop 

end 
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APPENDIX B. A NUMERICAL 
PROCEDURE FOR INTERPOLATION OF 
IRREGULARLY DISTRIBUTED DATA 

This appendix describes the procedure used to interpolate data known at locations 
that are irregularly distributed throughout a plane. Given the finite set of known data 
values (x„ y t , z,) the interpolating function f(x, y) satisfies z,- = /(x,, y,) for all 
known values. The interpolating function is used to calculate interpolated values at 
the desired (x, y) coordinates. The term irregularly distributed means that it is not 
possible to construct a quadrilateral grid using the known data coordinates as the 
vertices of the quadrilaterals. 

Bivariate interpolation procedures can be categorized as global methods or local 
methods. Global interpolation methods use all the known data to calculate the 
interpolated value at the desired coordinate. Global methods can become excessively 
cumbersome when the number of known data points is large and can lead to unwanted 
oscillations in the interpolating function. Local interpolation methods offer an 
improvement to the problems of global methods by only using a subset of the known 
data to calculate the interpolated value at the desired coordinate. This requires a 
rational procedure to determine which of the known data are used to calculate the 
interpolated value. Often, this is done by partitioning the irregularly distributed data 
into a triangular grid using the known data coordinates as the vertices of the triangles. 
The procedure of constructing a triangular grid from the irregularly distributed data 
is referred to as triangulation. 

There is available in the IMSL library of Fortran subroutines a local bivariate 
interpolation program for irregularly distributed data. This program was written 
by Akima and is described in References 68, 69. Akima’s program constructs a 
triangular grid and calculates local interpolating functions that possess continuous 
first derivatives. This procedure was not used for two reasons. First, because the 
interpolating functions have continuous first derivatives there will be some oscillation, 
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although much less severe than global methods. Second, and more importantly, the 
numerical procedure described in Appendix C requires an intermediate result of the 
interpolation procedure, the triangular grid. This information cannot be obtained from 
the IMSL subroutine and the Akima program is copyrighted and the source code is 
unavailable. For these two reasons a new interpolation program was written. 

The procedure described in this appendix to interpolate irregularly distributed data 
was be divided into two parts. The first part involves constructing the triangular grid 
from the irregularly distributed data. The second part is the actual interpolation. 
In practice, the triangulation requires more sophisticated programming while the 
interpolation requires more computational effort. Each part is described in a separate 
section. 


Triangulation 

The triangulation procedure used was based on the methods described by Akima 
[68] and Lawson [70]. Triangular grids constructed from irregularly distributed data 
are not unique. Therefore, an optimizing criterion was used to determine which 
triangles are constructed. The optimizing criterion was designed to avoid creating 
long and thin triangles whenever possible. 

The set of known data is arranged in ascending order according to Euclidian 
distance from an arbitrary reference point. The first data point was used here as the 
reference point. Lawson used the point having the smallest x coordinate and Akima 
used the midpoint of the two nearest points. The first two data points of the ordered 
data are always vertices of the first triangle constructed. The third vertex of the first 
triangle is found by examining the ordered data sequentially for the first noncolinear 
point. The data are then reordered by placing the point used as the third vertex of the 
first triangle in the third position of the ordered data and increasing by one the indices 
of all points from the previous location of the third point to the previous location of 
the point used as the third vertex. 
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Each triangle is recorded by the three indices associated with the points forming 
the vertices. All triangles are recorded in counter-clockwise order. Additionally, 
a record is made of the indices of the points which form the perimeter of the 
region of all constructed triangles. The perimeter record is also maintained in 
counter-clockwise order proceeding successively around the perimeter. The program 
proceeds sequentially through the ordered data from the fourth through the last point 
constructing all possible new triangles using the new point as one vertex and points on 
the perimeter as the other two vertices of each new triangle. Proceeding sequentially 
through the ordered data insures that new points will lie outside the perimeter of 
the triangularized region. Gearing all possible triangles using the new point and 
perimeter points insures the triangularized region will remain convex. 

For a new point, the test used to determine which points on the perimeter are used 
to construct new triangles can be visualized by drawing a line from the new point to 
each point on the perimeter. If this line cuts through the perimeter then that particular 
perimeter point is not used. The actual procedure used involves calculating the angles 
formed by the new point and successive pairs of adjacent points on the perimeter. 

After all new triangles have been constructed for a new point, each new triangle 
and the neighboring triangle which shares a common edge are examined according 
to the local optimizing criterion. If the criterion requires rearranging the new triangle 
and its neighbor then two new triangles are constructed and each are examined for 
rearrangement. This process continues until all new triangles have been examined 
and all satisfy the local optimizing criterion. 

The local optimizing criterion decides which of two possible triangles is con- 
structed from the quadrilateral formed by two triangles sharing a common edge. The 
criterion used here calculates the four interior angles of the quadrilateral. The opti- 
mized triangles are constructed by placing the common edge through the point on the 
quadrilateral with the largest interior angle. The criterion used by Akima [68] and 
Lawson [70] is referred to as the max-min angle criterion which creates the pair of 
triangles with the greatest least interior angle. 
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Linear Interpolation 

In this program the interpolated value is calculated at the desired point (x,y). 
This is done by first determining which triangle die desired point is a member of 
and then calculating the interpolated value for the point Ii the point lies outside the 
convex hull containing all the known points, and therefore is not a member of any 
triangle, an arbitrary value (here zero) is recorded for the interpolated value. 

A method suggested by Akima [68] was used to accelerate the procedure which 
determines of which triangle the desired point is a member. This method is most 
useful when interpolated values are desired at many locations. This method avoids 
testing the entire set of triangles by dividing the region in which the interpolation is 
desired into four quadrants. A list is prepared for each quadrant of the triangles that 
have a nonempty intersection with that quadrant For each point where interpolation 
is desired it is first determined to which quadrant the point is a member. Then, only 
the list of triangles for that quadrant is tested to determine of which triangle the point 
is a member. 

After the correct triangle is found for the desired point the linear interpolation is 
calculated. To calculate the interpolation in the ith triangle, a basis B, given by Equa- 
tion B.l, is first calculated from the three points (xs,i,y,-,i), (xi,2,yi,2), 3 ) 

which are the vertices of the triangle. The coordinates (a, /?) of the point (x, y) with 
respect to the basis B are then calculated using Equation B.2. This result is used to 
calculate the interpolation using Equation B.3. 


B 

= {(x 2 - £1,3/2 - yi ),(£3 - 

ys - yi)} 

(B.l) 

a 

r i-i 

£«',2 £»,1 £f ,3 £*, 1 

X - X t) l 

(B. 2 ) 

[P , 

[ yt, 2 - y», 1 y.,3 - y*,i . 

. y - y«, 1 . 



2 = Zi,i + a(zi,2 ~ Zi, 1) + $(zi t 3 

- *», l) 

(B. 3 ) 
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A Fortran Program For Triangulation 

c ___ | | I I I I — | -c 

c c 

c TRIANGLE 2.0 (08/07/89) c 

c c 

c Bruce Reichert c 

c Mail Stop 5-11 c 

c NASA Lewis Research Center c 

c 21000 Brookpark Road c 

c Cleveland, Ohio 44135 c 

c (216) 433-8397 c 

c c 

c I I -I Purpose | | I I -c 

c c 

c The purpose of TRIANGLE is to triangulate the x-y plane (divide c 

c the x-y plane into triangular cells) using the supplied set of c 

c irregulary distributed 2D data points as the vertices of the c 

c triangular cells. These triangles can be used later for analysis c 

c of the data (e.g. interpolation, graphical display, statistical c 

c analysis) . c 

c c 

c- I |- Input - Output I -I |-c 

c c 

c The input to TRIANGLE is a file containing the set of 2D data c 

c points to triangulate. The first record must be the number of c 

c points in the dataset* The following records must contain the c 

c the x and y values in the first and second fields. Only the first c 

c two fields of each record are read, subsequent fields are c 

c ignored. c 

c c 

c The first record of the output file is the number of triangles c 

c constructed. The following records each have three fields c 

c containing the three vertices of a triangle. The vertices are c 

c written not as x-y coordinates but as the integer number of the c 

c datapoint (the order of the data points occurance in the input c 

c file) . c 
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c c 

c | | — r* Variables I j |-c 

c c 

c Real c 

c c 

c xydata - The array of data points used to construct the triangles, c 

c xydata is arranged in order from least to greatest c 

c distance from the the first data point. c 

c distan - The array containing the distance from the first data c 

c point. This is used only, to arrange the order of xydata. c 

c • c 

c Integer c 

c c 

c index - The array that contains the original location of the data c 

c points in the input file. c 

c irow - The current data point in the array xydata. c 

c ldata - The . total number of data points in the array xydata. c 

c triang - The array containing the vertices of the triangles. c 

c mdata - The current number of triangles in the array triang. c 

c ntri - The number of new triangles constructed using the current c 

c data point. c 

c perim - The array containing the vertices of the triangles that c 

c lie on the perimeter of the triangulated domain. c 

c ndata - The current number of vertices in the array perim. c 

c ' c 

c Character c 

c c 

c filel - The name of the input file containing the set of data c 

c points. c 

c file2 - The name of the ouput file containing the vertices of the c 

c triangles. c 

c c 

c I I — • — I Method -| i — I | -c 

c c 

c To triangulate the x-y plane the data set is first ordered from c 

c least to greatest distance from the first point of the data set. c 

c Triangles are then constructed by proceeding sequentially through c 

c this ordered array. Ordering the data points from least to c 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

C- 


greatest distance insures that the new data point will lie outside c 
the domain of existing triangles. c 

c 

New triangles are constructed by determining what vertices of c 

existing triangles the new data point can "see", Only points on c 
the perimeter of the domain of the existing triangles need to be c 
considered. c 

c 

The newly constructed triangles are then inspected to determine c 

if they need rearranged. This is done by finding the existing c 

triangle that lies adjacent to the new triangle. The four c 

vertices of these two adjacent triangles form a quadralateral . c 

By constructing a diagonal line between opposite vertices of the c 

quadralateral two possible sets of triangles can be formed. c 

The triangles which are construced are those formed by placing c 

the diagonal through the largest interior angle of the c 

quadralateral. The vertices of triangles are always written in c 

counter-clockwise order. c 

c 


l-c 


common 

xydata. 

triang. 

perim 

common 

P it 

twopi, 

ntri 

common 

Idata, 

mdata, 

ndata 

common 

irow 




real xydata (10000,2) , distan ( 10000 ) 

integer triang (20000, 3 ) , perim ( 10000) , index ( 10000 ) 

character* 80 filel, file2 


pi = atan2(0.0, -1.0) 
twopi = 2 * pi 

write (*, ' (a)') ' Enter the cartesian data file name' 

read (* , ' (a) ' ) filel 

open (unit = 1, file = filel, status— 'old') 

write (*, '(a)') ' Enter the triangle data file name' 



192 


read (*, '(a)') file2 

open (unit = 2, file = file2, status = 'new') 
read (1, '{!)') ldata 

c Order the data from the file filel into the array xydata from 
c least to greatest distance from the first point in the file 
c filel. 

read (1, ' (2f ) ' ) x, y 


xydata (1, 1 ) = x 
xydata (1, 2) = y 
index (1 ) =1 

distan(l) = 0.0 


c Read the ith row of data from filel and calculate its distance 
c from the first point. 

do 10 i = 2, ldata 

read <1, '(2f)') x, y 

r = (x - xydata (1, 1)) ** 2 + <y - xydata (1, 2)) ** 2 

c Compare the distance of the point in the jth row of the array 
c xydata (which currently has values in rows 1 to i - 1) with the 
c distance of the point in the ith row of the file filel. If the 

c distance of the point in the ith row from filel is less than the 

c distance of the point in the jth row of the array xydata then 

c move all the values of the array xydata from rows j+1 to the end 

c (row i - 1) ahead one row to make room in the array xydata for 
c the new row. If the distance of the point in the jth row of the 
c array xydata is greater than or equal to the distance of all the 
c points in the rows of the array work than j will exit the 
c comparison loop with a value j = i and row shifting loop will 
c not execute. 

do 20 j = if i - 1 

if (r .It. distan(j)) goto 30 
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20 continue 

30 do 40 k = i f j + 1, -1 

xydata (k, 1) = xydata (k-1, 1) 
xydata (k, 2) = xydata (k-1, 2) 
index(k) - index(k-l) 

distan(k) = distan(k-l) 

40 continue 

c Now add the ith row of from file filel into the proper kth row 
c of the array xydata. 


xydata ( j, 1) = x 
xydata (j, 2) = y 
index(j) = i 
distan(j) = r 


10 continue 

irow = 2 

50 irow - irow + 1 

c Calculate the angle formed by the first two points seen from the 
c irow-th data point. Use this to determine if the three points can 
c form a triangle. 

theta = atan2 (xydata {2 , 1) - xydata (irow, 1), 

1 xydata(2, 2) - xydata(irow, 2)) 

2 - atan2 (xydata ( 1, 1) - xydata(irow, 1), 

3 xydata (1, 2) - xydata(irow, 2)) 

if (theta, le . 0.0) theta = theta + twopi 
if (theta .eq. twopi) goto 50 
if (theta, le. pi) then 

c Construct the first triangle and perimeter in the proper counter- 
c clockwise order. 
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triangd, 1) = 1 
triangd, 2) = 3 
triangd, 3) = 2 
perim(l) = 1 
perim(2) =3 
perim(3) = 2 

else 

triangd, 1) = 2 
triangd, 2) = 3 
triangd, 3) =1 
perim(l) = 2 

perim(2) = 3 

perim(3) = 1 

end if 

c Remove the irow-th values from the array xydata and index and move 
c rows 3 to irow - 1 ahead one position. Place the original irow-th 
c values in row 3. 

x = xydata (irow, 1) 
y = xydata (irow, 2) 
n = index(irow) 

do 60 j = irow, 4 , - 1 

xydata(j, 1) = xydata(j - 1, 1) 

xydata(j, 2) = xydata(j - 1, 2) 

index(j) = index(j - 1) 

60 continue 

xydata (3, 1) = x 
xydata (3 , 2Y = y 
index (3) = n 


mdata = 1 
naat a — 3 
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c Continue sequentially through the array xydata adding new triangles 
c and then rearranging them if necessary. 

do 70 irow = 4, ldata 
call newtri 
call arrang 

70 continue 

c Write the results to the output file. 

write (2, ' (i) ' ) mdata 
do 80 i = 1^ mdata 

write (2, '(3(i))') index (triang (i, l)) f index (triang (i, 2)) , 
1 index (triang (i, 3)) 

80 continue 

stop 

end 


subroutine newtri 


A 


c I I I Purpose | I I 

c 

c This subroutine determines what points on the perimeter the 
c new data point can "see" and then creates new triangles with 
c those points. After this is complete the perimeter array 
c is updated. 


c 

c I I Variables j I 

c 

c Real 
c 

c angle - An array that contains the angle between the new data 
c point and each point on the perimeter, 

c theta - The angle formed formed by two adjacent points on the 
c perimeter as seen from the new data point. 


l-c 

c 

c 

c 

c 

c 

c 

l-c 

c 

c 

c 

c 

c 

c 

c 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c- 


Integer 


c 

c 

c 


ifalse - The location preceeding the first occurance of the logical c 
value true in the array see. c 

itrue - The location preceeding the first occurance of the logical c 
value false in the array see. 

ntrue - The total number of occurances of the logical value true 
in the array see. 

work - An array used as a work space to update the array perim. 


Logical 

see - An array that contains a logical value of true for each 
point on the perimeter if it can be seen from the new 


data point. 

i i i i i 

i 

common 

i i 

xydata, 

i 

triang, 

i 

perim 

common 

Pi, 

twopi, 

ntri 

common 

ldata. 

mdata, 

ndata 

common 

irow 



real 

xydatad'0000, 

2 ) , angle (10000) 


integer 

triang (20000, 

3) , perim (10000) , 

work (1C 

logical 

see ( 10000 ) 




c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

I -c 


c Calculate the angle from the new data point to the points 
c on the perimeter 


do 10 i = 1, ndata 
j = perim ( i ) 

angle (i) = atan2 (xydata ( j, 1) - xydata {irow, 1) , 
1 xydata(j, 2) - xydata(irow, 2)) 

10 continue 


* 


k 


ntrue = 0 
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ntri = 0 

do 20 i = l f ndata - 1 

c Calculate the angle formed by two adjacent points on the 
c perimeter as seen from the new data point. Use this to 
c determine if the two points can be seen from the new data 
c point. 

theta ~ angle (i + 1) - angle (i) 
if (theta, le . 0.0) theta = theta + twopi 
if (theta, le . pi) then 


c Construct the new triangle. 


ntri 

mdata 

triang (mdata, 1 ) 
triang (mdata, 2) 
triang (mdata, 3) 


= ntri + 1 
= mdata + 1 
= perim(i) 

= irow 

= perimii + 1) 


ntrue 
see ( i) 


= ntrue + 1 
= .true. 


else 


see (i ) = . false . 

end if 
20 continue 

theta = angle (1) - angle (ndata) 
if (theta .le. 0.0) theta = theta + twopi 
if (theta . le . pi) then 


ntri 

mdata 


- ntri + 1 
= mdata + 1 
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triang (mdata, 
triang (mdata, 
triang (mdata, 


1) = perim(ndata) 

2) - irow 

3) = perim(l) 


ntrue 
see (ndata) 

else 


= ntrue + 1 
= .true. 


see (ndata) 


= .false. 


end if 

c Place the new perimeter points into the array work. 

if (see (1) ) then 

do 30 i = 2, ndata 

if ( .not .see (i) ) goto 40 
30 continue 

40 itrue = i - 1 

do 50 i = 1, ndata - ntrue 
work(i) = perim(i + itrue) 

50 continue 

if (itrue .It. ntrue) then 

work(i) = perim(i + itrue) 
else 

work (i) = perim(l ) 
end if 

i = i + 1 

else 

do 60 i = 2, ndata 
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if (see ( i) ) goto 70 
60 continue 

70 ifalse = i - 1 

do 80 i = 1, ndata - (ifalse + ntrue) 

work(i) = perim(i + (ifalse + ntrue)) 

80 continue 

do 90 i - ndata - (ifalse + ntrue) + 1, ndata - ntrue + 1 
work(i) ~ perim(i - (ndata - (ifalse + ntrue))) 

90 continue 

end if 

work ( i ) = irow 

ndata = ndata - ntrue + 2 

c Write the results of the new perimeter contained in the array 
c work into the array perim. 

do 100 i = 1, ndata 
perim(i) = work(i) 

100 continue 

return 

end 


subroutine arrang 

c I j j Purpose I J I I — c 

c c 

c This subroutine inspects the newly constructed triangles and c 

c rearranges them if necessary. If triangles are rearranged c 

c they are listed as new triangles and inspected for rearrangment . c 

c c 

c I j Variables I | |-c 
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c c 

c Real c 

c c 

c alphax - This is the angle from the xth vertex to adjacent vertex c 
c moving in the counter-clockwise direction c 

c betax - This is the interior angle of the the xth vertex of the c 

c quadralateral . c 

c c 

c Integer c 

c c 

c pointx - This is the xth point of the quadralateral . c 

c c 

C I : 1 — - I [ I — I — — I -c 

common xydata, triang, perim 

common pi, twopi, ntri 

common ldata, mdata, ndata 

common irow 

real xydata ( 10000, 2} 

integer triang (20000, 3), perim(lOOOO) 

integer pointl, point2, point3, point4 


c Select the new triangle. 

10 do 20 i = mdata - ntri + 1, mdata 

pointl = triang(i, 1) 
point2 = triang(i, 2) 
point3 = triang(i, 3) 

c Find the adjacent triangle. 

do 30 j = mdata - ntri, 1, -1 

do 40 k = 1, 3 

if (pointl .eq. triang(j, k)) then 
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do 50 1 - 1, 3 

if (point3 .eq. triang< j, 1) ) then 
c Form the quadralateral with the adjacent triangle. 

points = triang(j, 6 - k - 1) 

c Calculate the angles between the vertices of quadralateral. 



alphal = atan2 (xydata (point2, 

1) 

- xydata (pointl , 

i), 

1 

xydata (point2, 

2) 

- xydata (pointl, 

2)) 


alpha2 = atan2 (xydata (point3 , 

1) 

- xydata (point2 , 

1), 

1 

xydata (point3. 

2) 

- xydata (point2, 

2) ) 


alpha3 = atan2 (xydata (point 4 , 

1) 

- xydata (point3 , 

1), 

1 

xydata (points f 

2) 

- xydata (point3, 

2) ) 


alpha4 = atan2 (xydata (point 1 , 

1) 

- xydata (point4 , 

1), 

1 

xydata (pointl, 

2) 

- xydata (point4 , 

2)) 


c Calculate the interior angles of the quadralateral . 

betal - alphal - alphas + pi 

if (betal .It. 0.0) betal = betal + twopi 

if (betal . ge . twopi) betal = betal - twopi 

beta2 = alpha2 - alphal + pi 

if (beta2 .It. 0.0) beta2 = beta2 + twopi 

if (beta2 .ge. twopi) beta2 = beta2 - twopi 

beta3 = alpha3 - alpha2 + pi 

if (beta3 .It. 0.0) beta3 = beta3 + twopi 

if (beta3 .ge. twopi) beta3 = beta3 - twopi 

beta4 = alphas - alpha3 + pi 

if (beta4 .It. 0.0) beta4 « beta4 + twopi 
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if (beta4 . ge . twopi) beta4 = beta4 - twopi 

c Rearrange the triangles if necessary. 

if (max(beta2, beta4) .gt. max(betal, beta3)) then 
triang(j f 1) = triang(mdata - ntri, 1) 
triang(j, 2) = triang(mdata - ntri, 2) 
triang(j, 3) = triang(mdata - ntri, 3) 
triang(mdata - ntri, 1) = pointl 
triangtmdata - ntri, 2) = point2 
triang(mdata - ntri, 3) = point4 
triang(i, 1) = points 
triang(i, 2) = point2 
triang(i, 3) = point3 
ntri = ntri + 1 
goto 10 
end if 

goto 20 

end if 

50 continue 

end if 

40 continue 

30 continue 

20 continue 


return 

end 
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A Fortran Program for Linear Interpolation 


4 


4 


C I I — I I * I I -I -c 

c C 
c LINEAR 2.1 (02/19/90) c 
c C 
c Bruce Reichert c 
c Mail Stop 5-11 c 
c NASA Lewis Research Center c 
c 21000 Brookpark Road c 
c Cleveland, Ohio 44135 c 
c (216) 433-8397 c 
c c 

c I j I Purpose I I I l-c 

c c 
c The purpose of LINEAR is provide a 2D linear interpolation between c 
c irregulary distributed 2D data points. c 
c c 

c I I Input - Output I I l-c 

c c 
c The inputs to LINEAR are a file containing the set of 2D data c 
c points, the file produce by TRIANGLE containing the triangles c 
c constructed from the set of 2D data points, and the file c 
c containing the x and y coordinates of the points where the c 
c interpolation is desired. c 
c c 
c The output is the interpolated value for the given location. It c 
c is written in the output file as f(i) for the point x(i), y(i). c 
c c 

c I I Variables I I I -c 

c c 
c Real c 
c c 
c xydata - The array of 2D data points and their function value to c 
c interpolate between. c 
c x, y - Two arrays that contain the x and y coordinates of the c 
c points where interpolation is desired. c 
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c rx, ry - Arrays that define the x and y coordinates of the corners c 
c of the four regions. c 
c f - The interpolated value at the point x(i), y(i) . c 
c c 
c Integer c 
c c 
c ldata - The total number of data points in the array xydata. c 
c mdata - The total number of triangles in the array triang. c 
c ndata - An area containing the number of triangles in each region c 
c triang - The array containing the vertices of the triangles. c 
c nxy - The total number of x and y coordinates in the arrays c 
c x and y. c 
c nreg - The array that contains which region the triangles are c 
c members of. c 
c ntri - The number of the triangle that the point x(i), y{i) c 
c lies within . c 
c region - An array that contains which triangles are members of each c 
c region . c 
c c 
c Character c 
c c 
c filel - The name of the input file containing the set of 2D data c 
c points. c 
c file2 - The name of the input file containing the set of triangles c 
c constructed from the set of 2D data points. c 
c file3 - The name, of the input file containing the set of x c 
c coordinates where interpolation is desired. c 
c file4 - The name of the ouput file containing the interpolated c 
c values . c 
c c 

c 1 I I Method - | I I I -c 

c c 
c The interpolation domain is first partitioned into four rectangular c 
c regions. The list of triangles is searched to find which triangles c 
c the points defining the corners of the four regions are members of. c 
c Then it is determined to which region the points that define the c 
c triangles are members. This produces a list of triangles that c 
c are members of each region. For the given point x(i), y ( i ) where c 
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c interpolation is desired it is first determined in which region the c 
c point lies. Then the list of triangles that are members of that c 
c region are searched to find to which triangle the point is a member. c 
c A linear interpolation is then calculated using the three c 

c vertices of the correct triangle. If the point lies outside the c 
c data domain the interpolated value is recorded as 0. c 

c c 

c I I I I 1 I t-c 

common xydata, triang, region 

common ldata, mdata, ndata 

real xydata ( 10000, 3 ) 

integer triang (20000, 3 ) , region (20000, 5) , ndata (5) 
real x(1000), y(1000), rx(4, 2), ry<4,2) 

character*80 filel, file2, file3, file4, filed 
logical member 

write (*, '(a)')' Enter the irregularly distributed 

+ cartesian data file name' 
read ( * , ' (a) ' ) filel 

open (unit = 1, file = filel, status = 'old') 
read (1, ' (i) ' ) ldata 

write (*, '(a)')' Enter the triangle data file name' 
read ( * , ' (a) ' ) f ile2 

open (unit = 2, file = file2, status = 'old') 
read (2, ' (i) ' ) mdata 

write (*, '(a)')' Enter the grid data file name' 
read (*, ' (a) ' ) file3 

open (unit = 3, file = file3, status - 'old') 
read (3, ' (i) ' ) nxy 

write {*, ' (a)')' Enter interpolation data file name' 
read ( * , ' (a) ' ) filed 

open (unit = 4, file = file4, status = 'new') 
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do 10 i - 1, nxy 

read (3, ' (lx, 2 (fl0.5) ) ' ) x(i), y(i) 

10 continue 

do 20 i = 1, ldata 

read (1, ' (f 14 . 6, f 12. 6, f!3. 6> ' > xx, yy r ff 
xydata(i, 1) = xx 
xydatad, 2) = yy 
xydata (i, 3) = f f 
20 continue 

do 40 i = 1, mdata 

read <2, ' (3i) ' ) 11, 12, 13 
triang(i, 1) = 11 
triang(i, 2) = 12 
triang(i, 3) = 13 
region (i, 5) = i 
40 continue 

ndata (5) = mdata 

c Determine the min and max of the new grid. 

xmin - x { 1 ) 
xmax = x ( 1 ) 
ymin = y { 1 ) 
ymax = y ( 1 ) 

do 45 i = 2, nxy 

if (xmin . gt . x(i)) xmin = x(i) 

if (xmax .It. x(i)) xmax = x(i) 

if (ymin .gt. y ( i ) ) ymin = y(i) 

if (ymax .It. y ( i ) ) ymax = y(i) 

45 continue 

c Determine the boundaries for the four regions. 
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xmid = (xmin + xmax) / 2.0 
ymid = (ymin + ymax) / 2.0 


rx (1, 

1) 

- 

xmin 

ry(l, 

1) 

= 

ymin 

rx (1, 

2) 

= 

xmid 

ry (I, 

2) 

= 

ymid 

rx (2, 

1) 


xmid 

ry (2, 

1) 

= 

ymin 

rx (2, 

2) 


xmax 

ry (2, 

2) 

- 

ymid 

rx (3, 

1) 

= 

xmin 

ry (3, 

1) 

= 

ymid 

rx (3, 

2) 


xmid 

ry (3, 

2) 

= 

ymax 

rx (4, 

1) 

= 

xmid 

ry { 4 , 

1) 

= 

ymid 

rx (4 , 

2) 

= 

xmax 

ry (4, 

2) 

= 

ymax 


c Find which triangles the 16 points that define the 
c four regions belong to. 

do 50 i - 1, 4 

call which(rx(i, 1), ry(i, 1) , 5, ntri) 
if (ntri ,ne. 0) then 

call check (ntri, i, member) 
if (member .eq. .false.) then 
ndata(i) = ndata(i) + 1 
region (ndata (i ) , i) = ntri 
endif 
endif 

call which (rx(i, 2), ry(i, 1), 5, ntri) 
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if (ntri .ne. 0) then 

call check (ntri, i, member) 
if (member .eq. .false.) then 
ndata(i) = ndata(i) + 1 
region (ndata (i), i) - ntri 
endif 
endif 

call which(rx(i, 1), ry(i, 2), 5, ntri) 
if (ntri ,ne. 0) then 

call check (ntri, i, member) 
if (member .eq. .false.) then 
ndata (i) = ndata (i) + 1 
region (ndata (i) , i) = ntri 
endif 
endif 

call which(rx(i, 2), ry(i, 2), 5, ntri) 
if (ntri .ne. 0) then 

call check (ntri, i, member) 
if (member .eq. .false . ) then 
ndata (i) = ndata (i) + 1 
region (ndata (i) , i) = ntri 
endif 
endif 

50 continue 

c Find to which regions the three points that define each 
c triangle belong to. 

do 60 i = 1, ndata(5) 
do 70 j = 1, 3 
1 = triang (i, j) 
xtri = xydata (1, 1) 
ytri = xydata (1,2) 
do 80 k = 1, 4 

if (xtri . ge . rx (k, 1 ) ) then 
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if (xtri .le. rx(k, 2)) then 
if (ytri .ge. ry(k, 1)) then 
if (ytri .le. ry(k, 2) ) then 
call check (1, k, member) 

if (member ,eq. .false.) then 
ndata(k) = ndata(k) + 1 
region (ndata (k) , k) = i 
endif 
end if 
end if 
end if 
end if 

80 continue 

70 continue 

60 continue 


do 90 i = 1, nxy 
do 100 j = 1, 4 

c Find which region the point belongs. 

if (x(i) .ge, rx(j, 1)) then 
if (x (i) . le . rx ( j, 2) ) then 

if (y ( L) .ge. ry(j, 1) ) then 
if (y (i) .le. ry( j, 2) ) then 

c Find which triangle in this region the point belongs to 

call which (x (i) , y(i), j, ntri) 

c Calculate the linear interpolation. 

if (ntri .ne. 0) then 

call linint (x(i), y(i), f , ntri) 
else 

f = 0.0 


end if 
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end if 
end if 
end if 
end if 

100 continue 

write (5, ' <f ) ' ) f 

90 continue v 

stop 

end 

subroutine which (x, y, nreg, j) 

c This subroutine determines which triangle the point (x,y) lies within. 

common xydata, triang, region 

common ldata, mdata, ndata 

real xydata { 10000, 3 ) 

integer triang (20000, 3) , region (20000, 5) , ndata(5) 

pi .= atan2(0.0, -1.0) 

twopi = 2 * pi 

do 10 i = 1, ndata(nreg) 

j = region (i, nreg) 
xl = xydata (triang (j, 1) , 1) 
yl = xydata (triang (j, 1) , 2) 
x2 = xydata (triang (j, 2) , 1) 
y2 = xydata (triang ( j, 2) , 2) 
x3 = xydata (triang (j, 3) , 1) 
y3 - xydata (triang (j, 3) , 2) 

c Check to see if the point (x,y) coincides with a vertex. 


4 
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if ( (x.eq.xl ) .and. (y .eq.yl) ) go to 20 
if ( (x.eq.x2) .and. (y .eq.y2) ) go to 20 
if ( (x.eq.x3) .and. (y ,eq.y3) ) go to 20 

c Check to see if the point (x,y) lies within the interior. 

thetal = atan2(xl-x, yl-y) 
theta2 = atan2 (x2-x, y2-y) 
theta3 « atan2(x3-x, y3-y) 

phil2 = thetal-theta2 
phi23 = theta2-theta3 
phi31 = theta3“thetal 

if (phil2 . It .0 . 0) phi!2 = phil2 + twopi 

if <phi23 . It .0.0) phi23 = phi23 + twopi 

if (phi31 . It . 0 . 0) phi31 = phi31 + twopi 

if ( ( (phi 12 . le .pi ) . and . (phi23 . le .pi) ) .and . (phi31 . le .pi) ) goto 20 

10 continue 

j = 0 

20 return 

end 

subroutine linint (x, y , f , i) 

c This subroutine calculates a linear interpolation at the point x,y using 
c the vertices of the ith triangle. 

common xydata, triang, region 

common ldata, mdata, ndata 

real xydata (10000,3) 

integer triang (20000 , 3 ) , region (20000, 5) , ndata (5) 
real bvl(3), bv2(3) 
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xl 

yi 

fi 
x2 
y2 
f2 
x3 
y3 
f 3 


xydata (triang (i, 

i). 

1) 

xydata (triang (i. 

i). 

2) 

xydata (triang (i. 

i). 

3) 

xydata (triang (i. 

2), 

1) 

xydata (triang (i. 

2), 

2) 

xydata .(triang (i, 

2), 

3) 

xydata (triang (i, 

.3), 

1) 

xydata (triang (i. 

3), 

2) 

xydata (triang (i, 

3), 

3) 


c Create a basis for the plane. 


b v 1 ( 1 ) 
bvl (2) 
bvl (3) 
bv2 (1) 
bv2 (2) 
bv2 (3) 


(x2 - xl) 
{y2 - yl) 
<f 2 - fl) 
(x3 - xl) 
(y3 - yl) 
(f 3 - fl) 


c Calculate the coordinate of the point (x, y) with respect 
c to the new basis. 


det = bvl (1 ) 

* bv2 ( 2 ) - 

bv2 (1) 

* bvl (2) 



cl = (<x - xl) 

* bv2 (2 ) - 

(y - yl) 

* bv2 (1) ) 

/ 

det 

c2 = ( (y - yl) 

* bvl (1 ) - 

(x - xl) 

* bvl (2) ) 

/ 

det 


f = fl + cl * bvl (3) + c2 * bv2(3) 


return 

end 

subroutine check (ntri, nreg, member) 

c This subroutine- checks if triangle ntri is already 
c listed with the triangles belonging to region nreg. 

common xydata, triang, region 

common Idata, mdata, ndata 
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real xydata ( 10000, 3 ) 

integer triang (20000,3) , region (20000, 5) , ndata (5) 
logical member 

^ do 10 i = 1, ndata (nreg) 

* if (ntri .eq. region (i, nreg) ) then 

member = .true, 
return 
end if 

10 continue 

member = .false, 
return 

end 
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APPENDIX C. A NUMERICAL PROCEDURE 
FOR CALCULATING MOMENTS 
OF A BIVARIATE DISTRIBUTION 

This appendix describes the procedures used to calculate the moments of a 
bivariate distribution. In practice the distribution is not known, only values of the 
distribution at a finite number of points are known. In the first section the integrals 
used to determine the moments are examined by first considering the case where* 
the distribution is known explicitly and then considering the actual case where the 
distribution is only known at discrete locations. The second section includes the 
Fortran source code which implements the algorithms derived in the first section. 

The Calculation of Moments 

The moments of the distribution / which were calculated are defined by Equation 
C. 1. The integrals are over the domain Q which contains the support of the distribution 
/. The integral of the distribution over Q does not necessarily equal one and therefore 
appears in the denominator of Equation C.l. 

UfdA 

E[4>) = » <f> = x,y,x 2 ,xy,y 2 (C.l) 

From these moments the means, variances and covariance of the distribution are 
found using Equations C.2 through C.6 

H x = E[x\ (C.2) 


Vy = E[y] 


(C.3) 


a 2 x = E 


( x - n x y 


= E[x 2 J - E 2 [x] 


(C.4) 
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o] = E[(y - (<,) 2 ] = £[!, 2 ] - E 2 fa] (C.5) 

<T xy = E[(x - Hx)(y - fly)) = E[xy ] - E[x]E[y\ (C.6) 


A 




The distribution must be known to evaluate these moments. In practice, / is only 
known at a finite number of points, corresponding with the locations of experimental 
measurements. Because of the limited knowledge of the distribution, the calculations 
described are actually estimates of the moments. However, the calculated values 
approach the moments, with some restrictions on /, after sufficient refinement of the 
partition. 

For a partition {A t }” =1 of the domain fl which satisfies the conditions C.7, 
the moments may be calculated as shown in Equation C.8. The partition used 
in calculating the moments is the one determined by the triangulation procedure 
described in Appendix B. The elements of the partition are the sets of points 
A , which are the interior of the ith triangle whose vertices are the three points 
(^i,y«,i),(®»,2,y«,2),(a?»,3,y*,3)- Values of / are known only at the vertices. 

a - us, sns-fi,. \i) (c7) 


m= 


J 4>fdA 

Q 


f fdA 
n 


E I <t>fdA 

»=1A, 

tj fdA 

i'=l Ai 


(C.8) 
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To calculate the integrals used in evaluating the moments the change of variables 
shown in Equation C.9 was used. These change of variables have the beneficial effect 
of changing the integration from the elements of the partition A, to the set D shown 
in Figure C.1 


a 


x.,2 - x;,i 

- .*», 1 

-1 

ll 

A 


. Vi, 2 - Vi, 1 

Vi, 3 “ Vi, 1 . 


iH 

1 

1 


(C.9) 


y 

4 


(xi,3,yi,3) 


0 

4 


► a 


Figure C.l The change of variables from A, to D 



A reasonable approximation of f to use in evaluating the integrals is the linear 
interpolation described in Appendix B. Using the change of variables described in 
the preceding paragraph, the linear approximation of the distribution is determined 
by Equation C.10. 


/ = /»,! + <*(/», 2 - fi, l) + P(fi , 3 - fi, l) 


(C.10) 



217 


By using the change of variables it was possible to calculate the integrals exactly 
for the linear interpolation of /. Equation C.11 is a result of the change of variable 
theorem of integration of a function of several variables. 




dA = dxdy = 


d(a,l3) 


dad/3 


(C.11) 


The Jacobian was evaluated using Equation C.12, where Ai is the area of the element 
of the partition A,-. 


d(x,y) 


Xi, 2 - *»,l 

Xi, 3 1 

d(a,l3) 


Vi, 2 - Vi, 1 

Vi, 3 - Vi, 1 



£»,1 

T 

' 0 1 - 1 ' 




Xi, 2 


-1 0 1 


Vi, 2 


. *i',3 . 


_ 1 -1 0 . 


_2/»,3. 


= 2Ai 


(C.12) 


Since the points (xi ) i,y,- i i),(x,- ) 2,y» ) 2) 5 (®t,3»y» 1 3) and the values of the distribution at 
those points /,-,i,/;, 2 ,/j ,3 are constant, only products of a and /? up to second order 
remain within the integral. Thus, the integrals were evaluated in advance without the 
specific knowledge of the partition or the distribution. All that remains is algebraic 
operations involving the points (x,,i, y,,i), (x,,2, Vi, 2 ), (z»',3, y«',3) and the values of the 
distribution at those points /*,i , /*,2, /i,3- 

An example calculation of the integral used to evaluate the moment E[x ] is shown 
in Equation C.13. The results of calculations for the remaining integrals are given in 
Equations C.14 through C.18. The coefficients aj, b } and Cj^j are given in Table C.2 
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J xfdA = J (xj,i + a(x,,2 - xi,i) + £(*»,3 - x »,l)) 

d(x,y) 


A, 


(/«,! + a {fi,2 — /l) + ftifi, 3 — /»,l)) 


d(<*,ft) 


dad/3 



’®;,i 

T 

r 

'(1 — a — /?) 2 

a(l — a — ft) 

ft(\-a-py 


/*,1 

2A,- 

x i ,2 

/ 

a(l — a — ft) 

a 2 

aft 

dad/3 

ft, 2 


. x «.3. 

y 

A 

ft(l-a-ft) 

aft 

ft 2 


L/aJ 



T 

‘ 1/6 1/12 1/12' 


/•’, l 

Z .,2 


1/12 1/6 1/12 


fi,2 

_ X *,3. 


1/12 1/12 1/6 


. fi,2 . 


= Ai 


3 3 

= ^iXX 6 ^^* 

i=i Jfe=i 


/ 

A, 


fdA — A,- X! a jf*J 
i=i 


(C.13) 

(C.14) 


3 3 

I yfdA = A, - ^ ^ ^ bjkyjfk 
i i=i *=i 


(C.15) 


3 3 3 


/ «! «/ «/ 

x 2 fdA = A t - X X X c iXl x mfi 

A, tel*=W=l 


(C.16) 


3 3 3 


/ 4I 4# 4# 

y 2 fdA = Ai X X X c hk,iyj x kfl 

A, I=1*=1>=1 


(C.17) 


3 3 3 


/ xyfdA = Ai X X X c hk,iVj x kfl 

A. 1=1 k= i j=l 


(C.18) 
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Table C.l The values of the coefficients ay, 

fy.fc and c?,*,/ 

a i 

i = l 

j = 2 

J=3 


1/3 

1/3 

1/3 

b j,k 

& = 1 

jfc = 2 

k = 3 

j = 1 

1/6 

1/12 

1/12 

j = 2 

1/12 

1/6 

1/12 

J = 3 

1/12 

1/12 

1/6 

c i,k,i 

l = 1 

1 = 2 

1 = 3 

k = 1 

1/10 

1/30 

1/30 

k = 2 

1/30 

1/30 

1/60 

k = 3 

1/30 

1/60 

1/30 

c 2,k,l 

Z = 1 

1 = 2 

1 = 3 

II 

h- 1 

1/30 

1/30 

1/60 

k = 2 

1/30 

1/10 

1/30 

k = 3 

1/60 

1/30 

1/30 

cz,k,l 

Z = 1 

1 = 2 

Z = 3 

k = l 

1/30 

1/60 

1/30 

k = 2 

1/60 

1/30 

1/30 

k = 3 

1/30 

1/30 

1/10 

A Fortran Program For Two-Dimensional Moments Calculation 


c I I I 1 I I I -c 


C 


c 

C 

STATISTIC 1.1 (08/17/89) 

c 

C 


c 

C 

Bruce A. Reichert 

c 

C 

Mail Stop 5-11 

c 

c 

NASA Lewis Research Center 

c 

c 

21000 Brookpark Road 

c 

c 

Cleveland, Ohio 44111 

c 

c 

(216) 433-8397 

c 
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c c 

c j I I purpose | — I — I 1 -c 

c c 

c The purpose of STATISTIC is to calculate the means, variances, c 

c and covariance of a continuous bivariate distribution. The set c 

c of data are discrete values of the distribution which may be c 

c located irregularly throughout the plane. c 

c c 

c I I — Input - Output I 1 |-c 

c c 

c There are two files required as input to STATISTIC. The first c 

c is a file containing the set of data. The first record must be c 

c the number of points in the data set. The following records must c 

c contain the x and y values in the first and second fields and the c 

c function value in the third field. Any remaining fields are c 

c ignored. c 

c c 

c The second input file required is the output of the program c 

c TRIANGLE. The first record is the number of triangles constructed c 

c from the data set. The following records each have three c 

c fields containing the three vertices of a triangle. The vertices c 

c are written as the integer number of the datapoint (the order c 

c of the data points occurance in the input file) . c 

c c 

c The output of STATISTIC is a file containing the mean of x, mean c 

c of y, variance of x, variance of y and covariance. c 

c c 

c I I Variables I I |-c 

c c 

c Real c 

c c 


c xydata - The array of data points taken from the first input file. c 
c cX - Arrays that are used to calculate the integerals for each c 


c individual triangle. * c 
c area - The areal of the current triangle. c 
c el - The integral 'E [ 1] for the current triangle. c 
c ex - The integral E(X] for the current triangle. c 
c ey ^ - The integral E[Y] for the current triangle. c 


* 


* 
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c exx - The integral E[X*X] for the current triangle. c 

c exy - The integral E[X*Y] for the current triangle. c 

c eyy - The integral E [ Y*Y] for the current triangle. c 

c sumei - The sum of each integral E[l] for all triangles. c 

c sumex - The sum of each integral E[X] for all triangles. c 

c sumey - The sum of each integral E [Y] for all triangles. c 

♦ c sumexx - The sum of each integral E[X*X] for all triangles. c 

c sumexy - The sum of each integral E[X*Y] for all triangles. c 

c sumeyy - The sum of each integral E[Y*Y] for all triangles. c 

c xmean - The mean x value of the distribution. c 

c ymean - The mean y value of the distribution. c 

c xvar - The x variance of the distribution. c 

c xycov - The xy covariance of the distribution. c 

c yvar - The y variance of the distribution. c 

c c 

c Integer c 

c c 

c ldata - The total number of data points in the array xydata. c 

c mdata - The total number of triangles. c 

c point - The array containing the three points the current c 

c triangle . c 

c c 

c Character c 

c c 

c filel - The name of the input file containing the set of data c 

c points. c 

c file2 - The name of the input file containing the vertices of the c 

c triangles . c 

c file3 - The name of the output file containing the statistics of c 

c the data set . c 

c c 

c — | j | Method -| i i — |-c 

* C c 

c This program calculates the statistics for each individual triangle c 
r c and then sums each statistic for all triangles. For each triangle c 

c a linear approximation is assumed and the statistic is calculated c 

c exactly. c 

c c 
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c — I I 1 1 1 I l-c 

real xydata ( 10000, 3), x(3), y(3), f(3) 

real cl (3), c2<3,3), c3(3 f 3), c4(3 f 3 f 3) 

integer point (3) 

character*80 filel, file2, file3 

c Establish some parameters to save some writting. 

parameter ( pi = 1.0 / 3.0, p2 = 1.0 / 6.0 ) 

parameter { p3 = 1.0 / 10.0,.p4 = 1.0 / 12.0 ) 

parameter ( p5 = 1.0 / 30.0, p6 = 1.0 / 60.0 ) 

c Create the arrays used to calculate the integrals. 

data cl / pi, pi, pi / 

data c2 / 0.0, 0.5, -0.5, -0.5, 0.0, 0.5, 0.5, -0.5, 0.0 / 
data c3 / p2, p4, p4, p4, p2, p4, p4, p4, p2 / 

data c4 / p3, p5, p5, p5, p5, p6, p5, p6, p5, 

1 p5 , p5, p6, p5 , p3 , p5 , p6, p5, p5, 

2 p5 , p6, p5, p6, p5, p5, p5 7 p5, p3 / 

write (*, '(a)') ' Enter cartesian data file name' 

read ( * , ' (a) ' ) f ilel 

open (unit = 1, file = filel, status = 'old') 

write (*, '(a)') ' Enter triangle data file name' 

read ( * , ' (a) ' ) f ile2 

open (unit = 2, file = file2, status = 'old') 

write (*, '(a)') ' Enter statistic file name' 
read (*, ' (a) ' ) file3 

open (unit = 3, file = file3, status = 'new') 

c Load the data set from the first input file into the array xydata. 

read (1, ' (i) ' ) ldata 

do 10 i = 1, ldata 
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read(l, '(3f)') xydata(i, 1), xydata(i, 2), xydata(i, 3) 
10 continue 

c Proceed sequentially through the list of triangles calculating 
c the necessary integrals. 

* read (2, '(i) f ) mdata 

do 20 i = 1, mdata 

read (2, '(3i)') point (1), point (2), point (3) 

do 30 j = 1, 3 

x ( j ) = xydata (point (j) , 1) 
y(j) = xydata (point (j) , 2) 
f(j) = xydata (point (j) , 3) 

30 continue 

c Set the integrals for the new triangle to zero 

area = 0.0 
el = 0.0 
ex = 0.0 
ey =0.0 
exx = 0.0 
exy = 0.0 
eyy = 0.0 

c Calculate the integrals for one triangle. 

do 4 0 j = 1 , 3 

el = el + f ( j ) * cl(j) 

r do 50 k = 1, 3 

area = area + x(j) * y(k) * c2(j, k) 
ex = ex + x(j) * f (k ) * c3(j, k) 
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ey = ey + y(j) * f(k) * c3(j, k) 
do 60 1 = 1, 3 


exx = exx 

+ 

X ( j ) 

* x (k) 

* f (X) 

* c4(j. 

k, 1) 


exy = exy 

+ 

X(j) 

* y(k) 

* f(l) 

* C4(j, 

k, 1) 


eyy = eyy 

+ 

y(j) 

* y(k) 

* f (1) 

* C4(j, 

k, 1) 

4 


60 continue 


50 continue 


40 continue 


area = 

= abs 

(area) 

el = 

el * 

area 

ex 

ex * 

area 

ey - 

ey * 

area 

exx = 

exx 

* area 

exy = 

exy 

* area 

eyy = 

eyy 

* area 


c Sum the integrals for all triangles, 

sumel = sumel + el 

sumex = sumex + ex 

sumey = sumey + ey 

sumexx = sumexx + exx 
sumexy = sumexy + exy 
sumeyy = sumeyy + eyy 

20 continue 0 

c Normalize the resulting sum of integrals. ^ 

xmean = sumex / sumel 
ymean = sumey / sumel 
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xvar = sumexx / sumel - xmean * xmean 
xycov = sumexy / sumel - xmean * ymean 
yvar = sumeyy / sumel - ymean * ymean 

write (3,'<6f)') xmean, ymean, xvar, xycov, yvar, sumel 

stop 

end 
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